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Abstract 


A rigorous stability estimate for arbitrary order of accuracy of spatial central difference schemes 
for initial-boundary value problems of nonlinear symmetrizable systems of hyperbolic conservation 
laws was established recently by Olsson and Oliger (1994) and Olsson (1995), and was applied 
to the two-dimensional compressible Euler equations for a perfect gas by Gerritsen and Olsson 
(1996) and Gerritsen (1996). The basic building block in developing the stability estimate is 
a generalized energy approach based on a special splitting of the flux derivative via a convex 
entropy function and certain homogeneous properties. Due to some of the unique properties of the 
compressible Euler equations for a perfect gas, the splitting resulted in the sum of a conservative 
portion and a non-conservative portion of the flux derivative, hereafter referred to as the “Entropy 
Splitting”. There are several potential desirable attributes and side benefits of the entropy splitting 
for the compressible Euler equations that were not fully explored in Gerritsen and Olsson. The 
paper has several objectives. The first is to investigate the choice of the arbitrary parameter that 
determines the amount of splitting and its dependence on the type of physics of current interest to 
computational fluid dynamics. The second is to investigate in what manner the splitting affects 
the nonlinear stability of the central schemes for long time integrations of unsteady flows such as 
in nonlinear aeroacoustics and turbulence dynamics. If numerical dissipation indeed is needed to 
stabilize the central scheme, can the splitting help minimize the numerical dissipation compared to 
its un-split cousin? Extensive numerical study on the vortex preservation capability of the splitting 
in conjunction with central schemes for long time integrations will be presented. The third is 
to study the effect of the non-conservative proportion of splitting in obtaining the correct shock 
location for high speed complex shock-turbulence interactions. The fourth is to determine if this 
method can be extended to other physical equations of state and other evolutionary equation sets. 
If numerical dissipation is needed, the Yee, Sandham and Djomehri (1999) numerical dissipation 
is employed. The Yee et al. schemes fit in the Olsson and Oliger framework. 
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I. Introduction 


The construction of efficient high order low dissipation numerical methods for nonlinear 
conservation laws has been the subject of much research recently. For smooth flows, it is well 
known that the standard high order non-dissipative central schemes generate spurious noise leading 
to nonlinear instability, especially for long time integration applications such as in aeroacoustics, 
rotorcraft dynamics and turbulence physics. On the other hand, central schemes in conjunction 
with linear numerical dissipations are too diffusive for the physical problems in question. At the 
same time the majority of the available high order high-resolution shock-capturing schemes are too 
CPU intensive for practical computations. In spite of their high-resolution capability for rapidly 
evolving flows and short term time integrations, for long time integrations these schemes often 
exhibit undesirable amplitude errors for aeroacoustics problems. Current focus has been mainly 
on algorithmic issues in constructing highly accurate methods away from boundaries. Rigorous 
stability estimates for accurate and appropriate numerical boundary conditions associated with 
fourth- or higher-order methods are equally important, and have been the major stumbling block 
in the theoretical development of these schemes for nonlinear systems. Most of the existing theory 
for nonlinear conservation laws and their finite discretizations is concerned with the initial value 
problem (IVP). Standard practice in computational fluid dynamics (CFD) involving boundary 
conditions relies on guidelines from theory for linear stability analysis of initial boundary value 
problems (IB VPs) or IVP theories with the boundary conditions ignored. These linearized stability 
guidelines are only necessary but not sufficient for nonlinear stability. Spatial nonlinear stability of 
IB VPs goes hand-in-hand with the appropriate amount of nonlinear numerical dissipation required 
to stabilize the spatial scheme. The delicate balance of the numerical dissipation for stability 
without the expense of excessive smearing of the flow physics after long time integrations poses a 
severe challenge for unsteady flow simulations of this type. Employing a nonlinear stable form of 
the governing equations in conjunction with the appropriate nonlinear stable scheme for IB VPs is 
one way of minimizing the use of numerical dissipation. 

Until recently it was not known how to derive the proper numerical boundary conditions for 
a rigorous stability estimate for conventional spatial high order central difference schemes 
for nonlinear hyperbolic IB VPs. Advances by Kreiss and Scherer (1977), Strand (1994) and 
Olsson (1995a) led to the construction of high order boundary difference operators that enabled 
the design of stable high order central schemes for linear hyperbolic systems. The major tool 
used to overcome the stumbling block is a generalized energy method. The energy method for 
deriving stability estimates for hyperbolic IB VPs was first applied to the nonlinear scalar case by 
Gustafsson and Olsson (1994) for both higher-order central and compact symmetric schemes. It 
was then generalized and extended to nonlinear systems of symmetrizable hyperbolic conservation 
laws by Olsson and Oliger (1994) and Olsson (1995), and it was applied to the two-dimensional 
(2-D) compressible Euler equations for a perfect gas by Gerritsen (1996) and Gerritsen and Olsson 
(1996). With these recent developments, renewed interest has emerged in the use of spatial central 
schemes where efficiency, simplicity and non-dissipative properties are their trademarks. See for 
example, Yee et al. ( 1999) and Yee and Sandham (1998). 
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1.1. Relevance and Motivation 

The basic building block in establishing a stable energy estimate for high order spatial central 
schemes for nonlinear hyperbolic conservation laws consists of two parts. The first is a special 
transformation of the conservation laws to an appropriate form for the application of the continuous 
energy estimate for a stable IB VP of the governing equations. The second is a compatible numerical 
boundary difference operator for the application of the discrete analogue of the continuous energy 
estimate for a stable IBVP of the discretized counterparts. The special transformation relies on 
the symmetnzability of the systems of nonlinear hyperbolic conservation laws, the possession 
of a convex entropy function and a suitable splitting of the flux derivative vector with certain 
homogeneous properties. The compatible boundary difference operator has to satisfy the discrete 
analogue of the integration-by-parts procedure used in the continuous energy estimate. Olsson and 
Oliger (1994) utilized the result of Harten (1983a) on symmetric forms for systems of conservation 
laws as the backbone. Convexity of the flux functions is not required. These building blocks in 
tum allow one to use a modified form of the energy estimate (or generalized energy estimate) in 
deriving a compatible set of numerical boundary conditions that are stable for the higher-order 
central differencing schemes. 

Olsson proved that conservation is possible for second-order central schemes. To obtain a 
rigorous estimate for higher-order central schemes, one must apply the scheme to the split form of 
the flux derivative, written in non-conservative form, in terms of the transformed variables. For 
the Euler equations, one can further simplify the final split form of the flux derivative leading to 
a conservative portion and a non-conservative portion. The non-conservative portion as well as 
the conservative portion can be recast in terms of the conservative variables, making it tractable 
for practical applications in CFD. The proportion of the conservative to non-conservative parts 
is dictated by a parameter which falls into two wide ranges of intervals. The resulting splitting 
is hereafter referred to as the entropy splitting of the flux derivative or entropy splitting for 
ease of reference. Here, the entropy splitting should not be confused with the traditional flux 
vector splittings such as the Steger and Warming splitting (1981) or other vanants. The traditional 
flux vector splitting splits the flux function into different parts and most often into upwind and 
downwind portions. However, the entropy splitting splits the flux derivative using the properties 
of the chosen entropy function and the symmetrizability of the conservation laws without reference 

to any upwinding. 

Harten showed that the viscous terms of the compressible Navier-Stokes equations can also be 
symmetrized. In this case, only symmetry is needed in the derivation of the energy estimate. Due 
to the parabolic nature of the boundary conditions, the homogeneity properties are not required for 
the Navier-Stokes equations. For the numerical study involving the compressible Navier-Stokes 
equations in the present study, we apply the entropy splitting to the inviscid fluxes and the 
symmetric form of the viscous terms is not used. This is an attempt to examine if entropy splitting 
of the inviscid flux derivatives alone will provide side benefits over the un-split approach. 


Active research in the use of the symmetric form of the governing equations was carried out 
by Hughes. Franca and Mallet (1986) and related recent work. Hughes et al. utilized only 
the symmetric idea and employed the physical entropy as one of the entropy variables. Their 
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resulting inviscid flux vector and transformed state vector are not homogeneous in the entropy 
variables. Unlike the entropy splitting, their transformed equations are in purely non-conservative 
form. They have enjoyed improved results over the standard conservation law formulation. Their 
approach, however, does not allow a rigorous stability estimate for IB VPs for nonlinear hyperbolic 
conservation laws. In addition, due to their use of the purely non-conservative form, it is not 
certain that a correct shock speed can be obtained in general. 

The entropy splitting is not limited to smooth solutions. Olsson and Oliger (1994) also extended 
their result to weak solutions (problems containing discontinuities) that are obtained as pointwise 
limits of vanishing viscosity solutions. The entropy equality condition for the smooth solution 
case now becomes an entropy inequality condition. In addition, appropriate numerical dissipation 
is needed in conjunction with central schemes to pick out the physically relevant solutions if 
weak solutions are present. Gustafsson and Olsson (1995) proposed a scalar filter as numerical 
dissipation. Gerritsen and Olsson (1996) proposed the use of a slightly different nonlinear scalar 
filter in conjunction with wavelets for sharp shock capturing and shock detection. The recently 
developed high order low-dissipative shock capturing schemes using characteristic filters of Yee 
et al. (1999) fit in the entropy splitting framework. Instead of applying a scalar filter, they supply 
nonlinear filters based on, locally, the different wave characteristics of the convective flux. For 
complex shock waves, shear and turbulence interactions, one has better control of the amount 
of dissipation associated with each wave. For efficiency, Yee et al. proposed a combination 
of a narrow *rid stencil higher-order non-dissipative classical spatial differencing schemes and 
low order total variation diminishing (TVD), essentially non-oscillatory (ENO) or weighted ENO 
(WENO) dissipations as nonlinear characteristic filters with an artificial compression method 
(ACM) sensor. The ACM sensor is the same as that of Harten (1978) but applied in a slightly 

different context. 

The design principle of the Yee et al. schemes consists of two steps. The first step is the high 
order spatial and temporal base scheme. Unlike the hybrid schemes, the base scheme is always 
activated. The high order central and compact symmetric schemes are their primary choice for 
base schemes. In other words, the primary base schemes used by Yee et al. are exactly the central 
schemes for which Olsson and Oliger ( 1994) have provided rigorous stability estimates. The second 
step is the appropriate filter for stability and to capture shocks, shear-layers and fine scale flow 
structure. Many of the TVD, positive, ENO and WENO dissipations, after a minor modification 
by the ACM sensor, are suitable candidates as filters. The final grid stencil of these schemes 
is five in each spatial direction if second-order TVD schemes are used as filters and seven if 
second-order ENO schemes are used as filters for a fourth-order base scheme. There is only a 10% 
increase in operations count over standard second-order TVD schemes for 2-D direct numerical 
simulations (DNS). Based on their preliminary study for shock-turbulence flows, higher accuracy 
was achieved with fewer grid points when compared with that of standard second-order TVD, 
positive or ENO schemes. See Yee et al. ( 1999) or references cited therein. For the suppression of 
unphysical high frequency oscillations due to insufficient grid resolution and nonlinear instability 
away from discontinuities, higher-order spectral-like filters (Vichnevetsky (1974), Lele (1994), 
Alpert (1981). Visbal and Gaitonde (1998), Gaitonde and Visbal (1999)) might be needed where 
the value of the ACM sensor is extremely small. 
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1.2. Objectives 


Motivated by the aforementioned development of entropy splitting, Yee et al. (1999) proposed, 
as a followup work, to apply their schemes to the entropy splitting form of the tnviscid flux 
derivatives. This paper is a sequel to Yee et al. Besides investigating some of the fundamental 
issues described below, studies will be conducted to determine to what extent the entropy splitting 
form of the flux derivative can help in minimizing numerical dissipation, or equivalently, in 
improving nonlinear stability in conjunction with the Yee et al. (1999) schemes. Our main goal is 
to explore the possible side benefits of the entropy splitting without considering the accompanying 
stable numerical boundary difference operator as a complete package for stability requirements. 
This is accomplished by choosing numerical examples with periodic boundary conditions, or 
computational domains whose boundaries are far enough away so as not to affect the mainstream 
flow activities, and/or by using lower order non-characteristic boundary schemes. 


There are several potential desirable attributes of the entropy splitting for the compressible Euler 
equations that were not fully explored in Gerritsen and Olsson. First, in regions of smooth flows, 
additional numerical dissipation might not be required by the entropy splitting in conjunction 
with non-dissipative spatial central difference schemes. Second, the splitting appears to improve 
nonlinear stability over the un-split approach employing the same non-dissipative higher-order 
central schemes even for periodic boundary conditions. Third, the non-conservative portion of 
the flux derivative seems to have a small effect in obtaining correct shock speeds on the physical 
problems that Gerritsen and Olsson considered. Fourth, the entropy splitting in conjunction with 
higher-order central differencing could be a good candidate for nonlinear aeroacoustics, rotorcraft 
dynamics and turbulence computations where simplicity, high accuracy and low numerical 
dissipation are essential. But most of all, the splitting could possibly be extended to other physical 
equations of state and other evolutionary equation sets. The objective of this paper is many 
fold. The first is to investigate the choice of the arbitrary parameter that determines the amount 
of splitting (conservative and non-conservative proportions) and its dependence on the type of 
physics of current interest to CFD. The second is to investigate in what manner the splitting affects 
the nonlinear stability of the central schemes for long time integrations of unsteady flows such as 
in nonlinear aeroacoustics and turbulence dynamics. If numerical dissipation indeed is needed to 
stabilize the central scheme, can the splitting help minimize the use of numerical dissipation over 
its un-split cousin? Extensive numerical study on the vortex preservation capability of the splitting 
in conjunction with the central schemes for long time integrations will be presented in section 
4.1. At present, existing finite discretizations (applied to the un-split approach) that are capable 
of preserving long distance vortex convection are ineffective and CPU intensive. Extensive grid 
refinement and grid adaptation, and an unusually demanding time step reduction are necessary. 
The third is to study the effect of the non-conservative proportion of the splitting in obtaining 
the correct shock location for high speed complex shock-turbulence interactions. The fourth is to 
determine if this method can be extended to other physical equations of state and evolutionary 
equation sets. 

Aside from stability considerations, as explained in Harten (1983a), another potential desirable 
attribute in the use of the symmetric form of the governing equations is for the computation of the 
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steady-state solution of the conservation laws. In solving the steady nonlinear conservation laws, 
the symmetry of the matrix coefficients could possibly enhance the structure of iterative matrices in 
direct Newton-iteration methods. For time-marching to steady states or time accurate subiteration 
procedures using implicit methods, the symmetric form in conjunction with the entropy splitting 
might result in an improved convergence rate over the un-split approach. This will be a subject of 
future research. 


Outline: 

Section II reviews the entropy splitting and the numerical schemes for the 2-D compressible 
Euler equations for a perfect gas. The choice of the entropy splitting parameter is discussed in 
Section 2.2. Section III describes the extension of the entropy splitting to other physical equations 
of states and evolutionary equation sets. Section IV illustrates the performance of the entropy 
splitting for a variety of unsteady flows and compares the results with those obtained using the 
un-split conservative approach (i.e., without splitting the inviscid flux derivatives). The study 
concentrates only on perfect gases. The comparison uses the same spatial and time discretizations 
for both the split form and un-split form of the inviscid flux derivatives. If numerical dissipation 
is required, the Yee et al. (1999) filters are used as the numerical dissipation for nonlinear stability 
and/or to insure physically relevant numerical solutions and to suppress spurious oscillations across 
discontinuities. If spurious high frequency oscillations are present, the use of spectral-like filters 
(Vichnevetsky (1974), Lele (1994), Alpert (1981), Visbal and Gaitonde (1998), Gaitonde and 
Visbal (1999)) in conjunction with the Yee et al. filters might be needed at the locations where the 
value of the ACM sensor is very small. See Section 2.6 for a discussion. 

In this paper, unless indicated, Euler or Navier-Stokes equations pertain to compressible fluids. 
High order central difference schemes refer to fourth or higher-order spatial central difference 
schemes (compact or non-compact methods) without numerical dissipation added. Compatible 
time discretizations (in terms of stability and accuracy) should be used, but these are not the subject 
of this paper. The terms “split” and “un-split” mean the application of the same discretization to 
the “entropy splitting of the inviscid flux derivative” and the standard “inviscid flux derivative 
in terms of the conservative variables without splitting”. 


II. Entropy Splitting for a Perfect Gas 

This section reviews the basic building blocks for the entropy splitting for the 2-D compressible 
Euler equations for a perfect gas. The mathematical theory is quite involved. Interested readers are 
referred to references cited. The Yee et al. (1999) numerical methods used in conjunction with the 
entropy splitting are also summarized. 

2.1. Summary of Entropy Splitting for a Perfect Gas 


In vector notation the 2-D compressible time-dependent Euler Equations in conservation form 
for an equilibrium real gas can be written, in Cartesian coordinates, as 


7 


Ut + F, + Gy = 0 , 


(2.1.1a) 


where U t = F a = and G y = ff and the U , F, G, are vectors given by 
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The dependent variable is the vector of conservative variables, and ( p , u, v,p) is the vector of 
primitive variables. Here p is the density, u and v are the velocity components, pu and pv are the 
x- and y-components of the momentum per unit volume, p is the pressure, e - p[e + ( u + v )/2] 
is the total energy per unit volume, and e is the specific internal energy. 


For a thermally perfect gas, the equation of state is 

p = pRT , (2-1-2) 

where R is the specific gas constant, and T is the temperature with e = e(T). For constant specific 
heats (calorically perfect gas) 


e = c„T, (2-1.3) 

where c v is the specific heat at constant volume. 

The eigenvalues associated with the flux Jacobian matrices of F and G are (u,u,u ± c) and 
(v,v,v ± c), where c is the sound speed. The two u,u and v,v characteristics are linearly 
degenerate. Hereafter, we refer to the fields associated with the u ± c and v ± c charactensttcs as 
the nonlinear fields and the fields associated with the u,u and v,v characteristics as the linear 

fields. 

Here we outline the entropy splitting for a perfect gas in 2-D Cartesian coordinates. Formulas for 
the corresponding 3-D case can be found in Appendix B. Gerritsen and Olsson ( 1996) extended the 
summation-by-parts idea of Strand (1994), and the entropy splitting of Olsson (1995) and Olsson 
and Oliger ( 1994) to the 2-D Euler equations for an ideal gas in conjunction with high order central 
schemes. The first step in deriving the entropy splitting for the compressible Euler equations for a 
perfect gas is to introduce a symmetry transformation from the vector of conservative variables U to 
a new vector of symmetry variables W, referred to as the ‘ ‘entropy variables . The transformation 
is chosen so that the flux Jacobian matrix with respect to W is symmetric, and the Jacobian matrix 
of the transformation is symmetric and positive definite. A family of symmetry transformations, 
based on a scalar convex function tj, referred to as an “entropy function”, derived for the Euler 
equations for a perfect gas by Harten (1983a), was employed by Gerritsen and Olsson. It has the 

form 
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7 = P((S). (2-1-4) 

The function £( 5 ) is an arbitrary but differentiable function of a dimensionless physical entropy 

5 = log {pp~ y ) (2.1.5) 

where 5 has been non-dimensionalize by c„. The entropy variables W arejhen given by W = §£. 
The entropy function rj is not to be confused with the “physical entropy’’ 5 or the entropy variables 
vector W. The next step is to restrict the transformations to those that allow a special splitting of 
the flux derivatives. This requires that both the flux and conservative variables are homogeneous 
functions of the new variables. Under these conditions one can rigorously establish a bound on 
the rate of growth of the energy norm in terms of the absolute eigenvalues corresponding to the 
in-coming characteristic variables at the boundary of the domain. The relation between entropy 
functions and symmetrizable variables can be found in Mock (1980). 

In other words, one introduces an entropy variable transformation W = W(U) such that 
F w = f£ and Gw = jw are symmetric, and U w = jw is symmetric and positive definite. 
The entropy variable W is chosen such that F(U{W)), G{U(W)) and U(W) are homogeneous 
functions of W of degree /3, i.e., there is a constant /3 such that for all r 


U{tW) = r^U(W), (2-1-6) 

F(tW) = t 0 F(W). (2-1-7) 

The homogeneity property implies that 

F w W = PF{U(W)) (2-1.8) 

U W W = PU. (2-1.9) 


Then the splitting of F t results in 


F ‘ = JTI (FwW) - + J7 i F wW ‘ ~ TTi f - + 13 + 1 


F, + ——F w W„ 0 / -1. (2.1.10) 


A similar splitting can be written for G y and U t . 

1 

For a perfect gas, the required entropy function is obtained by letting £(5) = Ke ^+^ , where 
K and a are constants. The corresponding W can be written as 
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w = [wi to 2 to, 10 *] T = y [c + 7 TTP ~P U H 


(2.1.11a) 


and the upper triangular part of the symmetric matrix Uw is 


Uw = 



a/?u apv 

a/ni 2 — p apuv 
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f p(u 2 +v 2 )-^p 
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(2.1.11b) 


The constants a and b are a = (1 - a - j)/a and b = t/(t - !)• Here > P* and P are re,ated 
through 


p* = xe r=^T5 = (2.1.11 c) 

with v = < 0. In the authors’ opinion, the simplest choice is to set K = (3. The parameter 0 

A. f) 

is given by 


0 = 


a_f 7 
1-7’ 


(2.1. lid) 


Using (7 1 lb), (2.1.2), (2.1.3), (2.1.11a) and (2.1.1 Id), we can show (see Section 3.1) that U, F 
and G are homogeneous functions of W of degree 0. The positive definite condition on U w (see 
Section 3.1) restricts a to the two ranges a > 0 or a < - 7 - 


The flux vectors, expressed in the W variables are given by 

fwm) = £(-“> -=?(»> + 4 tp')] t . (21Ue) 

G(U(W)) = p[-w, ^ =Sr+p' + (2.1. uf) 

The upper triangular parts of the symmetric matrices F(U(W))w and G(U{W))w . expressed in 
the U variables, are given by 



apv} — p 
u(apu 2 - 3p) 


apuv 

v(apu 2 - p) 
ti{apv 2 - p) 


«[f p(u 2 +v 2 )-bp] 

- b + cpu 2 + fp(u 2 + u 2 )u 2 - |p( u2 + v 2 ) 
u»[cp+ | p(u 3 + t> 2 )] 

u[6c^- + cp(u 2 + v 2 ) + fp(u 2 + U 2 ) 2 ] 

(2.1.118) 
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apv 2 — p v[\p(u 2 + u 2 ) - bp] 

u(apv 2 — p) uv[cp + §p(u 2 + v 2 )] 

v{apv 2 - 3p) -b*f + cpv 2 - |p(u 2 + u 2 ) + f p(u 2 + v 2 )v 2 

v[&c£ + cp(u 2 + u 2 ) 4- fp(u 2 + v 2 ) 2 ] 

(2.1. llh) 


2.2. Choice of the Entropy Splitting Parameter 

From the structure of (2. 1.10), the entropy splitting divides the flux derivative into a conservative 
and a non-conservative part. The ratio between the conservative and non-conservative parts depends 
on the choice of the parameter 0. Both Harten ( 1983a) and Gerritsen and Olsson (1996) introduced 
the parameter a. In the authors’ opinion, the introduction of a is not necessary. However, to adhere 
to the discussion when referencing their work, we retain the use of a in the perfect gas case. The 
convexity condition on the entropy function r\ restricts the value a to two possible ranges; namely, 
a > Oora < -7 (or equivalently,/? < or/? > 0). Although Gerritsen and Olsson considered 
the a < -7 range which Harten (1983a) overlooked, they set a = 1 - 27 (J3 = 1) in conjunction 
with high order central differencing schemes in all of their numerical examples. This choice of 
a corresponds to the splitting of the flux derivative into equal conservative and non-conservative 
proportions. They did not give any guidelines or examples of the effect of the choice of a on the 
quality of the numerical solutions for different flow physics. In addition, all of their examples deal 
with at most simple shock waves, if present. Wavelets are used as shock detectors and to guide 
the grid adaptation. Due to the type of problem they addressed and the dense clustering of the grid 
points near the shocks using very small time steps, it is not certain that correct shock speeds were 
really obtained with a reasonably practical time step and grid distribution. It is the purpose of this 
section to discuss the choice of the a parameter value. We discuss a > 0 and a < -7 separately. 

The a >0 (or/3 < j-^)case: 

This is the only case that Harten considered. This corresponds to a negative /? which results in 
a conservative proportion fr c = ^ > 1 and a non-conservative proportion fr nc = ^ < 0 . 
As a -» 0 + , fr c -V 7“ and fr ne -> (1 - 7) + . Here, the superscripts “+” and indicate 
the values approach the limit from above and below respectively. Thus, it appears that a > 0 
is “non-standard” or “nonphysical” in CFD practice in the sense that a larger than 100% of 
the conservative proportion and a negative non-conservative proportion is used. As a -* 0 + , the 
proportion becomes extremely unphysical. As a — » 00, fr e — > 1 + and fr ne — > 0~. Therefore as 
a — ► 00, the proportion becomes more physical. 

The a < -7 (or 0 > 0) case: 

The a < —7 case corresponds to a positive (3 and consequently, fr c < 1 and fr nc < 1. 
This range of a conforms more closely to standard CFD practices. We have the following five 
situations. 


G w 


' apv apuv 

v(apu 2 — p) 


P* 


where c = (1 — 27)7(7 — 1 ). 
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i) As a — > - 7 , 0 -» 0 + , fr c -* 0 + and fr nc -» 1 

ii) For a = 1 — 27 , 0 = 1, fr c - 1/2 and fr ne = 1/2 

iii) For 1 — 27 < a < — 7 , fr c <1/2 and fr nc > 1/2 

iv) Fora < 1 - 27 , fr c > 1/2 and fr nc <1/2 

v) As a -+ - 00 , fr e -* 1" and fr nc -* 0 + 

Section IV gives a parameter study of a for three different types of flow physics. 

2.3. Numerical Methods 


The spatial discretizations for weak solutions proposed in Yee et al. (1999) consist of two 
parts, namely, a base scheme and a filter. When numerical dissipations or filters are not used, 
the scheme consists of only the base scheme. This section discusses the base scheme and Section 
2.4 discusses the form of the filter (numerical dissipation) for complex shock waves, shear and 
turbulence interactions. Section 2.6 discusses the blending of the Yee et al. (1999) filters with 
other filters for the suppression of spurious high frequency oscillations. 

2.3.1. Spatial Base Schemes for the Convection Terms 

Denote Fj, h as the discrete approximation of the convection flux F at (j A®, k Ay), where A* 
and A y are the grid spacing in the x- and y-directions and j and k are the corresponding spatial 
indices. Possible non-dissipative high order base schemes for F a (similarly for G y ) can be of the 
following two types. 


Central Diffenncings : (fourth and sixth-order) 

F m as 12 {^ i+2,k ~ ^Fj+i.k 4- - Fj- 2,k^ , 

F, as _ + 45F;+i,fc — 45F ; _i,i, -I- 9Fj-2,h — Fj-3,uj- 


(2.3.1) 

(2.3.2) 


Compact Ccntval Differcncings'. (fourth and sixth-order) 


F m 



where for a fourth-order approximation 

(A a F)j,h = - (Fj+i,h + 4 Fj t j, + 


(2.3.3a) 

(2.3.3b) 

(2.3.3c) 
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and for a sixth-order approximation 

(A.F)i,u = + 3 Fj, u + (2.3.3d) 

(B,F ) m , = 1 (F i+hi + 28 Fj+i.k - 28 f+,.» - (2.3.3.) 


2.3.2. Spatial Schemes for Viscous Terms 

For simplicity let V aa be a viscous term in one dimension. The possible high order base schemes 
for V aa can be 

Central Differencing ». (fourth and sixth-order) 


12Aaj 2 


(Vj+2 ~ 16 V i+l + 30V) - 16^-! + V)_ 2 ) , (2.3.4) 

V„ « — i-y f 2V) + 3 - 27V)+ 2 + 270V)+i - 490V) + 270V)_x - 27V)_ 2 + 2V)_ , ) . (2.3.5) 
180Aas \ ' 


Compact Central Differencingy. (fourth and sixth-order) 


V aa *-^(c: l D a v) , (2.3.6a) 

where for a fourth-order approximation 

(C.V), = i (v j+1 + 10 V, + V,-.). (2.3.6b) 

(D.V)j=V j+l - Wj + Vj-u (2.3.6c) 

and for a sixth-order approximation 

{C.V), = V)+i + a 0 V) + V)-i, (2.3.6d) 

(■ D a V)j = bo ( Vj+i - 2V) + V)-x) + y ( v )+ 2 - 2V) + V)_ 2 ), 

V / 4 V / (2.3.6e) 

a 0 = 5.5, (2.3.6f) 

feo = 4(ao — 1 )/3, (2.3.6g) 

cq = (10 — a 0 )/3. (2.3.6h) 
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2.4. Filters 

In this section we first review the procedure for applying the characteristic filter to multistage 
Runge-Kutta type and linear multistep method (LMM) types of time discretizations (Yee et ah 
(1999)). Examples of explicit LMMS are forward Euler and Adams-Bashforth. Examples of 
implicit LMMs are backward Euler, trapezoidal rule, and three-point backward differentiation. 
The one-leg formulation of the LMMs of Dahlquist (1979) is also applicable. We then discuss 
forms of the characteristic filter. 

2.4.1. Procedure to Apply the Filter Step 


If a multistage time discretization such as the Runge-Kutta method is desired, the spatial 
differencing base scheme discussed in the previous section is applied at every stage of the Runge- 
Kutta method. If viscous terms are present, we use the same order and type of base scheme for the 
viscous terms as for the convection terms. 

There are two methods for applying the characteristic filter. Method 1 is to apply the filter 
at every stage of the Runge-Kutta step. Method 2 is to apply the filter at the end of the full 
Runge-Kutta step. For inviscid and strong shock interactions, method 1 might be more stable. 

If one desires a time discretization that belongs to the class of LMMs, then the filter can be 
applied as a numerical dissipation vector in conjunction with the base scheme. The filter in this 
case is evaluated at U n for explicit LMMs. For implicit LMMs additional similar filters evaluated 
at the n + 1 time level are involved. Alternatively, method 2 can be applied to LMMs as well. 
In this case, we apply the filter after the completion of the implicit time step. One can minimize 
flux evaluations by using the one-leg formulation of the LMMs of Dahlquist (1979). The only 
non-dissipative (in time) second-order, two-time level one-leg method is the mid-point implicit 
method. Note that the noniterative linearized form of the midpoint implicit formula reduces to the 
regular noniterative linearized trapezoidal formula. 

For time marching to steady states using implicit LMMs, certain flow physics only requires an 
explicit dissipation term. Also, the implicit operator can be different from the explicit operator. 
See Yee (1985, 1986, 1989), Yee et al. (1990) for some efficient conservative linearized implicit 

forms. 


2.4.2. Nonlinear Characteristic Filters 

There are many possible candidates for the filter operator in conjunction with high order base 
schemes. Here, we propose using filter operators that have similar width of grid stencils as the base 
scheme for efficiency and ease of numerical boundary treatment. Higher than third-order filter 
operators are of course applicable, but they are more CPU intensive and require special treatment 
near boundary points for stability and accuracy. 

We use nonlinear dissipation terms in conjunction with the Harten ACM sensor applied to each 
characteristic wave as the filter vector. In essence, the nonlinear dissipation terms act as second or 
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third-order ACM-like operators instead of Harten’s first-order ACM (Harten (1978)). The sensor 
is used to signal the amount of nonlinear dissipation to be added to the high order non-dissipative 
scheme, one wave at a time. Thus, the current approach is also different in spirit from using ACM 
to sharpen the contact discontinuities in the original Harten (1983b) second-order TVD scheme. 
Let the filter vector in the ^-direction be of the form 


F* is the modified form of the nonlinear dissipation portion of the standard numerical flux. 

c dF - 

For characteristic based methods, the quantity R j+ 1 is the right eigenvector matrix of w using, 
for example, Roe’s approximate average state. Note that the eigenvector R jJr i should not be 

confused with the R in (2.1.2). We cast the G) h+ , in the same manner. The elements of 


denoted by 



( 2 . 4 . 2 ) 


(f> 1 . , in (2.4.2) are the elements of $ j+ 1 - the dissipative vector of the high resolution schemes 
resulting from using a TVD, ENO or WENO scheme. Hereafter, we refer (2.4.2) as the ACM filter. 

The function kQ 1 . + j is the key mechanism for achieving high accuracy of the fine scale flow 
structure as well as shock waves in a stable manner. In other words, the elements of are the 
same as the nonlinear dissipation term of the TVD, ENO or WENO scheme with the exception 
of premultiplying by k6 1 . + , . The parameter zc is problem dependent. For smooth flows, k is 
used to improve nonlinear stability and can be very small. Different physical problems require a 
different value of «c because of the large variation in flow properties. The * value may vary from 
one characteristic wave to another, and from one region of the flow field to another region with 
different flow structure. The range of k for our present numerical experiments is 0.0 < * £ • 
The function 9 l . x is the Harten ACM sensor. For a general 2m + 1 points base scheme, Harten 

recommended 




1 

3 



— m+H 



( 2 . 4 . 3 ) 





la' 








( 2 . 4 . 4 ) 


The a' , are elements of R.^tiUj+i ~Uj)- 

j ■'’ r 5 

Instead of varying k for the particular flow problem, one can vary p. For larger p, less numerical 
dissipation is added. Note that by varying p > 1 in (2.4.4), one can essentially increase the order 
of accuracy of the dissipation term. The order of the dissipation depends on the value of p (Bjorn 
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Sjogreen, private communication, 1998). One can switch fromp 
at smooth regions. For all of the numerical examples, we use p - 




= max 




: 1 near shock locations to p > 1 
1 and 

(2.4.5) 


The shock-turbulence interaction problems appear to favor this form of 9 \ + , . 


Formulae for <f> 1 . ± are well known and can be found in the literature. In most of the numerical 
computations in Section IV, we use the Harten and Yee (1985) second-order upwind TVD 
numerical dissipation. Computations using the symmetric TVD dissipation (Yee, 1985) will also 
be presented. See Yee et al. (1999) for details and for a discussion of other possible filters. 


2.5. Computer Implementation 

To avoid some conditional statements in the actual computer code and to promote vectorization, 
several of the functions inside the filter with the potential of dividing by zero are modified. See 
Yee et al. (1999) for details. In particular, the sensor (2.4.4), with p = 1, is replaced by 

al _ IK+i/zl ~ ' a i-i/ 2 : (2.4.6) 

In all of the computations, we take e = 10~ 7 . (Actually, e should have the same dimension as 
a j+i/2^ 

2.6. Blending ACM Filters with Other Filters 


The ACM filter (2.4.2) might not be sufficient for (a) time-marching to steady state and (b) 
spurious high frequency oscillations due to insufficient grid resolution and nonlinear instability 
away from discontinuities, especially for turbulent and large-eddy simulations. This Section 
discusses the blending of other filters with ACM filters. 

(a) Time-Marching to Steady State: 

For time-marching to steady state one usually needs to add fourth-order linear dissipation to a 
second-order spatial differencing scheme (Beam and Warming (1976)). For the present schemes 
using characteristic filters, in addition to the ACM filter, one might need to add a sixth-order 
linear dissipation to a fourth-order spatial base scheme and an eighth-order linear dissipation to a 
sixth-order spatial base scheme in regions away from shocks for stability and convergence. Let L d 
be such an additional filter operator. Take the case of a Runge-Kutta time discretization. There 
are again two ways of incorporating the L d operator. One way is to incorporate the L d operator 
at every stage of the Runge-Kutta method. Another way is to include the L d operator after the 
completion of the Runge-Kutta full step. The best way of applying the L d operator is most likely 
problem dependent and time integrator dependent. For LMM type of time integrators, L d is used 
in conjunction with the ACM filter step as an additional dissipation. 
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To minimize the amount of dissipation due to L d in the vicinity of shock, waves, there should be 
a switching mechanism n d (different from /c in (2.4.2)) to turn off the L d operator in the vicinity 
of shock waves. The L d operator can be applied to the conservative, primitive or characteristic 
variables. The simplest form is to apply L d to the conservative variables. Alternatively, since 
all of the work in computing the average states and the characteristic variables is done for the 
shock-capturing filter operator, one can apply the L d operator to the characteristic variables. In 
this case, the switching mechanism k d can be a vector so that it is more in tune with the ACM 
shock detector using the approximate Riemann solver. For example, one can set k = 0 for the 
linearly degenerate fields and blend a small amount of * d to remove spurious noise generated by 
the lack of ACM filters. This blending of the ACM filter with the L d operator can be applied to 
time-accurate computations as well. 

(b) Suppression of Spurious High Frequency Oscillations: 

The ACM filters might not be able to remove spurious high frequency oscillations effectively 
unless sufficient grid points are used. For the suppression of unphysical high frequency oscillations 
due to insufficient grid resolution and nonlinear instability away from discontinuities, higher-order 
spectral-like filters (Vichnevetsky (1974), Lele (1994), Alpert (1981), Visbal and Gaitonde (1998), 
Gaitonde and Visbal (1999)) might be needed at the locations where the value of the ACM 
sensor is very small or zero. If spectral-like filters are needed, a proper blending of ACM filters 
with spectral-like filters should be applied. In this case, we can use the same procedures as the 
time-marching to the steady state except the L d operator should be replaced with the spectral-like 
filters (for compact central schemes). 

III. Extension to Other Equations of State and Equation Sets 

In this section we discuss the extension of entropy splitting to more general cases. The method 
originally was developed for the 2-D Euler equations in Cartesian coordinates for a perfect gas. 
We show here how it can be extended to flow of a gas that is only thermally perfect. For 
maximum generality, the analysis is presented for arbitrary three-dimensional, time varying grids. 
For compactness, we employ the vector approach of Vinokur (1989). Here the word vector 
refers to a physical vector such as velocity or momentum, as distinguished from an algebraic 
vector representing a set of variables. For completeness, the Roe Riemann solver for a thermally 
perfect gas is included. This is motivated by the fact that if the characteristic type of nonlinear 
dissipation or filter is desired, Roe’s Riemann solver is normally employed. It is noted that 
both Steger- Warming flux-vector splitting and Roe’s approximate Riemann solver have an exact 
extension for this case. For the readers unfamiliar with the vector approach of Vinokur (1989), the 
results for the 2-D and 3-D Euler equations in Cartesian coordinates are given in Appendices A 
and B. The caloric equation for an ideal diatomic gas is given in Appendix C. We also examine 
the possible extensions to a nonequilibrium mixture of species, magnetohydrodynamics, and the 
artificial compressibility method for incompressible flow. 


3.1. Entropy Splitting for a Thermally Perfect Gas 


In this subsection we derive entropy splitting for a gas that is only thermally perfect, with 
the internal energy being an arbitrary function of temperature. This law is valid for a dilute gas 
consisting of a single chemical species, and is also a very good approximation for air below the 
temperature when oxygen starts to dissociate (approximately 2000° K). 


The following development describes the derivation leading to a final form of the entropy 
splitting of the flux derivative for a thermally perfect gas. It has the same form as the perfect gas 
case, but, the corresponding ranges of the /3 parameter are different, and W, U w and F w have 
different expressions. In fact, the derivation of entropy splitting for a perfect gas has to be modified. 
Certain parameters that are constant for a perfect gas are no longer constant for a thermally perfect 
gas. In particular, Harten and Gerritsen normalized their entropy 5 by c„, which is now a variable 
for a thermally perfect gas. We therefore normalize our S by the gas constant R. This results in 
differences in our results from theirs involving the quantity y - 1. The positive definite condition 
(or equivalently, the convexity condition) on Uw again restricts to be in two possible ranges. As 
mentioned previously, Harten (1983a) overlooked one of the more physical ranges, and needlessly 
introduced a parameter a in his solution. Such a parameter has no analogue for the more general 
thermally perfect gas case, and only serves to complicate the derivation. Harten also introduced 
an arbitrary constant K , which he then set equal to a particular value in order to simplify the 
final expressions. We chose a particular value from the beginning, and avoided introducing an 
un-needed constant. 

We repeat the equation of state for a thermally perfect gas (2.1.2) 

P = P RT, (3.1.1) 


where p, p, T, and R are the pressure, density, temperature, and gas constant, respectively. The 
entropy 5 and internal energy e are then related to p and T through the first and second laws of 
thermodynamics by 

dS=i-^, (3.1.2) 

T P 


where we have introduced the normalized variables 

and T = RT. 




(3.1.3) 


Equation (3.1.1) can then be rewritten as 


p = P T. 


(3.1.4) 


While 5 is dimensionless, T has the same dimension as e. It follows from^(3.1.2) that e - e(T) 
only. All real gases satisfy the conditions e > 0 and c > 0, where k = dtj dT . 

Equation (3.1.2) can be integrated to obtain 
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where . 

/(f) = exp(-y (3.1.6) 

The arbitrary constant in the integral of (3.1.6) can be absorbed in the definition of 5. 

Since conservation laws are expressed in terms of of conserved quantities per unit volume, it 
is convenient to introduce the internal energy per unit volume c = pe. If u is the fluid velocity 
vector, then the set of conservative variables U can be represented by the vector 

U = [p m e] T , (3.1.7) 


where m = pu is the momentum vector per unit volume, and e — e + ^-pu • u is the total energy 
per unit volume. Note the algebraic vector U has three elements, of which the second element is 
the physical vector m. The temperature T(U) is obtained by solving the equation 



1mm 

2 ~?~' 


(3.1.8) 


Equation (3.1.8) has a unique solution since e > 0. Let n be the unit normal vector in a positive 
direction to a cell surface in a finite-volume grid, or a coordinate surface in a finite-difference 
grid. If v n is the normal component of the velocity of a time-varying surface, and u„ = u • n, one 
can define the normal relative velocity component u' = u n — v n . The set of inviscid normal flux 
components F n is given by the vector 


F n = [pu' mu'+pn eu' + pu n ] T . (3.1.9) 

Following Harten (1983a), we obtain the transformed variable W from 

(3.1.10) 

8U 

where the function rj(U) is given by 

T] = pZ(S). (3.1.11) 

Here again the second element of the algebraic vector F n is a physical vector quantity. £(S) is 
an arbitrary function of S , and 5 can be expressed as a function of U using (3.1.5), (3.1.6) and 
(3.1.8). Using (3.1.2) and (3.1.8), we can express the differential of (3.1.11) as 

drj = -|[e - 2e -p(l +0 )}dp - m • dm + pdej, (3.1.12) 


m 

where ( = d£/dS and 



(3.1.13) 
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Note that fl is in 
multiplication as 


general a function of S. We can rewrite (3.1.12) in the form of a matrix 
© 


dr,= - [e — 2? — p(l + /?) 

P 


-m 


P 


dp 

dm , 
de 


(3.1.14) 


where the vector dot product is implied in multiplying the second element of the row vector by the 
second element of the column vector. In the rest of the paper, we will use the convention that in 
forming the product of two matrix elements, a dot product is implied if each element is either a 
physical vector or a tensor. From (3.1.10) we can express drj as drj = W T dU. We thus obtain 


W — [w w w\ T = -[e — 27— p(lT/3) — m p] • (3.1.15) 

In order to investigate the homogeneity condition, it is useful to introduce the function 

p(W) = 2-^ = — e(f ) - f (1 + 0). (3-1-16) 

w 2(tF) 

Note that p involves only thermodynamic variables, and is also a homogeneous function of W of 
degree 0. Since e is in general a nonlinear function of f, the homogeneity condition can only be 
satisfied if f is a homogeneous function of W of degree 0. In view of (3.1.5) and (3.1.16) this can 

only be accomplished if ^ 

/3 (5) = constant. (3.1.17) 

We can now solve for £(5) from (3.1.13). The sign of £ determines the positive definite condition 
on Uw'< the scale of £ does not affect any numerical calculations. Anticipating the positive definite 
condition which will be derived below, we find that the simplest solution of (3.1. 13) is 


t =/3e-^, (3-1-18) 

which then gives 

£ = -e-*/ 0 . (3.1.19a) 

Substituting (3.1.5), this can be expressed as 

( = (3.1.19b) 

The expression for p(W) obtained using (3.1.4), (3.1.7), (3.1.15) and (3.1.19b) is 

(3-1-20) 

f(TV 

Therefore p is a homogeneous function of W of degree /3. From (3.1.15) and the definition of m 
it follows that u, u n and u' are all homogeneous functions of W of degree 0. We conclude from 
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definitions (3.1.7) and (3.1.9) that the vectors U and F n are both homogeneous functions of W of 
degree /3. 

While we cannot obtain an explicit formula for U(W ), we can derive an explicit expression for 
the symmetric matrix Uw as functions ol U . The upper triangular part can be written as 


Uw = 7 


' ap apu ae + bp 

apu u - pi u[ae + (b — l)p] 

^(i +j3). 


(3.1.21a) 


where 


«(?,/ 3 ) = 


1, 


i+1+0 

t(?,0) = (l+ l 9)a+g = - + 


— e 


(3.1.21b) 


(3.1.21c) 


and I is the identity tensor. 

Let A = Uw and u = qc, where c is a unit vector in the direction of u and q is the magnitude 
of u. The positive definite condition for A is determined by the quadratic form X T AX , where the 
vector X has the general form X = [ z\ *je *3 j r , and the dot product is implied in a product 
involving the unit vector e. Since e • ee • e = e*I*e = e*e = l, the quadratic form can be 
written in terms of ordinary scalar quantities as X ,T A' X 1 , where X’ = [ x\ *3 x s ] > and the 
upper triangular part of A 1 becomes 


( 


ap 


apq 

,2 


ae + bp 

apq * - p q[ae + (b - l)p] 

Therefore the positive definite condition for Uw is obtained by calculating the signs of 
principal minors of (3.1.22a). Since the elementary operation of subtracting a multiple 
from another leaves a determinant unchanged, we can reduce (3.1.22a) by a series of 
operations to the matrix 

” ap apq ae + bp 

0 - p -pq 

0 0 


1 

i 


pa 


(3.1.22a) 


the leading 
of one row 
elementary 


(3.1.22b) 


This matrix is positive definite if the leading principal minors are positive. Since p > 0 andp > 0, 
we obtain the conditions 

a/£ >0, (3.1.23a) 

—all 2 > 0, (3.1.23b) 


and 

b(3/i 3 > 0. (3.1.23c) 

From (3.1.23a) and (3.1.23b) it follows that l < 0, a condition already satisfied by (3.1.19a). It 
then follows from (3. 1.23c) that 6/3 < 0, which can be reduced to 

1 -1 

p > 177 ' 


(3.1.24) 
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Condition (3.1.24) is satisfied if (3 > 0 or 0 < -(1 + c). It is easy to show that a value of 0 in 
these two ranges satisfies a < 0, as required by (3.1.23a). Since e > 0, the maximum value of e 

occurs at T maa . Therefore, for (3 < 0, 

j3 < — [l + e(f m ..)]. (3.1.25) 

A sufficiency condition, independent of the flow problem, is obtained by replacing e(T moa ) by 
e(oo). 

We can also derive an explicit expression for the symmetric flux Jacobian with respect to the 
transformed variables, {F n ) w , as functions of U. The upper triangular part can be written as 

' u'ap u'apn — pn u'(ae + bp) - u n p 

(F n )w = ~ u'{apuu - pi) - p(un + nu) o 2J , (3.1.26a) 

( a i3 


where 


and 


a 2 , = {u'[ae + (b- l)p] - u„p}u - ^(e + p)n, 

— u 1 — - +pf — - u> u^) + — (1 + 0) -2 u„-(e + p). 

. p V p J p J p 


(3.1.26b) 

(3.1.26c) 


3.2. Roe Riemann Solver for a Thermally Perfect Gas 

The extension of Roe’s approximate Riemann solver to a thermally perfect gas has been given 
by Abgrall (1990) and also Spekreijse and Hagmeijer (1990). (They actually considered the more 
general case of a mixture of thermally perfect gases, valid for non-equilibrium flow.) We present 
the results here for arbitrary three-dimensional grids, using our compact vector notation. The 
Riemann solver is based on properties of the ordinary flux Jacobian matrix A = . From (3. 1 .7) 

and (3.1.9) it follows that we need the pressure differential 


dp 

= x dp + k de, 

(3.2.1) 

where - 

and x(T) = T — ne. 


<T) = j 

(3.2.2) 

The matrix A can then be written as 



-v n 

n 0 

(3.2.3) 

A = K\i\ - u n u 

un - xnu + u'l /cn , 

_(i U-H)u n 

H n — KU„U U 1 + KU n 


where K\ = • u + x- H = ^ + | u * 

u is the total enthalpy per unit mass, and h = 

e + T is 


the specific enthalpy. The three distinct eigenvalues of A are 

A 1 = u\ A 3 = u + c, and A 3 = u - c, 


(3.2.4) 


22 


where the speed of sound c is given by 

c 2 — x + “h- 


(3.2.5) 


The multiple eigenvalue A 1 is associated with those conservative variables whose flux is 
purely convective. These are p and the tangential component of m. In order to construct the 
corresponding linearly independent eigenvectors associated with the multiple eigenvalue, we span 
the plane normal to n by an arbitrary set of two basis vectors b< and the set of reciprocal 
basis vectors b^, satisfying b; • b^ = gj, where 6j is the Kronecker delta. It follows that 
b . . n _ b i . n _ o. The right eigenvector matrix R can then be written as 


1 


R = 


u 

K 2 


0 1 
cbi u + cn 
cbj • u H + cu n 


1 

u — cn , 
H - cu n 


(3.2.6) 


where K 2 = |u • u - x/ K - 

Among the various approximate Riemann solvers, the most common one uses the Roe average 
because of its simplicity and its ability to satisfy the jump conditions across discontinuities exactly. 
In those solvers based on local linearization, the flux at a point separating two states U and U is 
based on the eigenvalues and eigenvectors of some average A. The optimum choice for A is one 

satisfying _ . , 

AF n = A &U, (3-2-7) 


where A(-) = (■)* — (-) 1 '. This choice of A captures discontinuities exactly. One way of obtaining 
A is to seek an average state U, which is a function of U L and U R , such that 


A = A{U). 


(3.2.8) 


Such a state is known as a Roe-averaged state. One can easily show that 


u = au L + (1 - a)u H 

(3.2.9) 

and _ „ 

H = ctH L +{ 1 - a)H R , 

(3.2.10) 

where ^ 

Cl — — — 

(3.2.11) 

1 + V> H /P C 

From the definition of H one then obtains 

h = H - \ u ■ u. 
2 

(3.2.12) 


The discrete form of (3.2.1) yields the relation 

X A p + k Ae = A p. 


(3.2.13) 
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The Roe-averaged sound speed is given by (3.2.5) as 


c 2 = x + 


(3.2.14) 


Equation (3.2.13) provides only one relation to determine x and «. Since x and * ^ functions of 
f only, the simplest assumption is that x and * depend only on f fl and f L . Eliminating p, using 
(3.1.4), we rewrite (3.2.13) as 

x(p R _ /) + n( p R e R - pV) = p R T R - p L T L . (3.2.15) 

Equating the coefficients of p R and p L on the two sides of (3.2.15) we obtain the relations 


K = 


AT 

a7 


(3.2.16a) 


and 


X = 


g.R'pIt _ i j' R 

Ac ' 


(3.2.16b) 


Equations (3.2. 16a, b) are replaced by (3.2.2) when AE — * 0. An important quantity in the 
approximate Riemann solver is the column vector R~ l AU . Its components are the jumps in the 
characteristic variables. It is given by 


where 


R~ 1 AU 


A p — Ap/c 1 

pb J • Au/c 

| ( Ap/c 2 + pn • Au/c) 
| (Ap/c 2 - pn - Au/c) 



3.3. Entropy Splitting for Other Equation Sets 


(3.2.17) 


(3.2.18) 


In this subsection we examine the possibility of applying entropy splitting to other equation 
sets. We first consider non-equilibrium flow, consisting of a mixture of different species, each 
obeying a thermally perfect gas law. The motivation is again the fact that this case can also be 
treated exactly with Steger- Warming flux-vector splitting and Roe’s approximate Riemann solver. 
We next consider the equations of magnetohydrodynamics, since there is much current interest in 
their solution. Finally, we investigate the artificial compressibility method applied to the solution 
of the incompressible equations. 

3.3.1. Non-equilibrium Flow 

In non-equilibrium flow, we consider a mixture of species, each obeying a thermally perfect 
gas law. The conservation law now takes the form 


U t + V • F = S, 


(3.3.1) 
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where S is a vector consisting of source terms for each species. The equation of state for species : 

15 p’ = pWT, (3-3.2) 

where p‘, p‘, and R' are the pressure, density, and gas constant for species i, respectively, and T 
is the temperature of the mixture, assumed to be in thermal equilibrium. If there are ns species in 
the mixture, the index i takes on values from 1 to ns. The entropy 5* and internal energy «' are 


then related to p‘ and T by 

dS R'T p { ’ 

(3.3.3) 

where we have introduced the normalized variable 


II 

(3.3.4) 

Note that we have not introduced a normalized T, since the R l are different for each species. It 
again follows from (3.3.3) that = e^T) only. All real species satisfy the conditions e‘ > 0 and 
e* > 0, where c* = dc'/dT. Equation (3.3.3) can be integrated to obtain 

i ti — 5* 

p f — e , 

(3.3.5) 

where 

f'(T) — ex p( — J ^y dT )- 

(3.3.6) 

The arbitrary constant in the integral of (3.3.6) can be absorbed in 

the definition of 5*. 

We can now define the density of the mixture, p, as 


II 

-M 
•* . 

(3.3.7) 

the pressure p as 

p = '£p‘ = p rt , 

i 

(3.3.8) 

where __ 

p R = Y j P' R ‘. 

i 

(3.3.9) 

and the entropy per unit volume, pS, as 


,s = £ys‘. 

(3.3.10a) 


l 

Using (3.3.4), the last equation can be written as 

pRS = y>'R'S\ 

i 


(3.3.10b) 
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where 

S - (3.3.11) 

R 

Note that since R is no longer a constant, S is not proportional to 5. Finally, we define the internal 
energy of the mixture per unit volume as 

e = pe=^pV. (3.3.12) 


The set of conservative variables U can be represented by the vector 

U = [% m e] r , (3.3.13a) 

where the vector % is defined as 

% = (p 1 p 3 ... P nt ). (3.3.13b) 

The temperature T(U) is obtained by solving implicitly the equation 

X>V(T) = «- (3.3.14) 

i 

where p is given by (3.3.7). Equation (3.3.14) has a unique solution since p l > 0 and > 0 . The 
set of inviscid normal flux components F n is given by the vector 

F n = ['Ru' mu'+pn eu' + pu n ] T . (3.3.15) 

The procedure to obtain the transformed variable W follows that of subsection 3.1. Equa- 
tions (3.1.10)-(3.1.13), and (3. 1 . 17)-(3. 1 . 19a) in that subsection are still valid. With the aid of 
(3.3.3), (3.3.7), (3.3.9), (3.3.10b), and (3.3.14) we obtain 

W = [W w Wf, (3.3.16a) 


where 

W = (ui 1 w 2 ... w nt ) , (3.3.16b) 


w i = — 


RT 


[-!TT(1 + S - 5*) -S3RT - e‘ + * = l,...,nj, (3.3.16c) 


w = — 


pRT ’ 


(3.3. 16d) 
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and 



(3.3.16c) 


It is again useful to define a set of functions of W that are homogeneous of degree 0, involving 
only thermodynamic variables, as 


uHW) = = - = -£‘T(1 + 5 - 5‘) - 0RT - c*' i = 1, (3.3.17) 

w 2(^) 3 


Equations (3.3. 16e) and (3. 1 .19a) can be combined to yield 

(-Wf(RTf=e- 3 . (3.3.18) 


In order to prove homogeneity, we let 


P i = a*(-w) 


fi 


(3.3.19) 


Substituting (3.3.5), (3.3.18) and (3.3.19) into (3.3.17), we obtain 


+ W + 1 - 0H*T) + ln(«7‘) = 0. 


(3.3.20) 


Combining (3.3.4), (3.3.9), and (3.3.10b) to eliminate p' , yields 


R = 


Li ai ' 


Substituting (3.3.5), (3.3.18), and (3.3.19) into (3.3.10b), we obtain 


0]n(RT) = 


Y,i aiRi 


(3.3.21) 


(3.3.22) 


Equation (3.3.20), (3.3.21), and (3.3.22) comprise a set of coupled non-linear equations for R , T, 
and a* as functions of p'. Since the /i‘ are homogeneous functions of W of degree 0, it follows 
that R, T, and a { are all homogeneous functions of W of degree 0. It is then easy to show that U 
and F n are homogeneous functions of W of degree 0. 

In order to obtain an expression for Uw, we can combine the differentials of (3.3.5), (3.3. 17) and 
(3.3.18) to express dp * as a linear combination of dw l , dw, dw, dR, and dT. From the differentials 
of (3.3.4), (3.3.5), (3.3. 10b), and (3.3.18) we find that dT equals a linear combination of dw, dR, 
and all the dp'. By eliminating dp from the differentials of (3.3.4) and (3.3.5) we see that dR is also 
equal to a linear combination of all the dp \ Therefore, obtaining Uw requires inverting a dense 
linear system. It would thus be difficult to establish the positive definite condition. Therefore the 
extension of the method to non-equilibrium flow is not practically feasible. If the homogeneity 
condition is not required, then one can use symmetry variables based on the physical entropy, as 
was shown by Chalot et al. ( 1990). . 
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3.3.2. Magnetohydrodynamics 

The set of conservative variables U for magnetohydrodynamics is represented by the vector 

U = [p m e B] r , (3.3.23a) 

where B is the magnetic field vector, and 

e = pe{T) + l B. (3.3.23b) 

2 p l 

The set of inviscid normal flux components F n is given by the vector 

pu' 

pu'u + n(p + |B • B) — 2? n B 
eu' + (p + |B • B)u n — B n ( u • B 
u n B - B n u 

where B n = B • n. 

Equation (3.3.23a) shows that B is an element of U, while (3.3.23b) shows that ^B • B is a 
term in e, which is also an element of U . It follows that both B and e cannot have the same degree 
of homogeneity, and therefore U cannot be a homogeneous function. This result is not surprising, 
since the magnetic field vector is not a true conservative variable, and is not expected to behave 
the same way as the physical conservative variables. 

3.3.3. Artificial Compressibility Method for Incompressible Flow 

For the artificial compressibility method for incompressible flows, the set of conservative 
variables U for this case is given by 

U = [p u] T , (3.3.25) 

while the set of inviscid normal flux components F n is given by 

F n = [cru 1 u'u-t-pn] T , (3.3.26) 

where <r is the artificial compressibility. Since the second element of U is u, while one of the 
terms in the second element of F n is uu • n it follows that U and F n cannot have the same 
degree of homogeneity. This is also not surprising, since p is not a true conservative variable. 
In the compressible case, the velocity has homogeneity of degree 0, and it is the density that is 
homogeneous of degree j3. When the density is no longer a variable, the homogeneity property 
disappears. 

IV. Numerical Examples 

The numerical experiments will be limited to a perfect gas. Three test cases are considered. 
The first is inviscid and the last two are compressible mixing layer computations. These test cases 



(3.3.24) 


28 


were also considered in Yee et al. (1999). The three test cases are: (1) a horizontally converting 
vortex, (2) a vortex pairing in a time-developing mixing layer with shock waves formed around the 
vortices and (3) a shock wave impinging on a spatially evolving mixing layer where the evolving 
vortices must pass through a shock wave, which in turn is deformed by the vortex passage. For 
the two mixing layer computations, the study will be limited to the choice of the $ that determines 
the amount of splitting in obtaining the same shock location as the un-split approach. For the 
Navier-Stokes computations involving entropy splitting, the splitting is applied to the invisc.d flux 
terms, and the symmetric form of the viscous flux is not used (see Section I) 


In all of the computations the classical fourth-order Runge-Kutta time discretization is employe 1 
For the purposes of this paper we concentrate on the non-compact central schemes (-3.1) and 
n 3 2) with the same order of accuracy and type of base scheme for the convection and viscous 
terms (if viscosities are present). Compact symmetric schemes (2.3.3) are also applicable^ 
but require nearly twice the CPU time that the non-compact central schemes require for 2-D 
compressible mixing layer computations and will not be addressed in this paper. 


If numerical dissipation is added, the filters (2.4.2), (2.4.5) and (2.4.6) are used at the end of 
the full Runge-Kutta time step. Roe’s (1981) average states are used in (2.4.1). For most ot the 
computations, the Harten and Yee (see Yee and Harten (1985a) and Yee (1989)) second-order 
upwind TVD dissipation for *J +1 in (2.4.2) is used. These will be notated as ACM with the 

following numbers indicating the order of the base scheme for the convection and viscous terms. For 
example, ACM44 means the use of fourth-order central as the base scheme for both the convection 
and viscous terms. In order not to introduce additional notation, inviscid flow simulations are 
designated by the same notation, with the viscous terms not activated. Computations using 
symmetric TVD dissipation (Yee, 1985, limiter (2.7b)) are indicated by adding the letter S ^ , as in 
ACM44S. Computations using entropy splitting are indicated by adding the letters ENT at t e 
end as in ACM44-ENT. To examine the performance of the entropy splitting schemes where shock 
waves are absent, the computations also employ only the non-dissipative central base schemes 
(without the ACM filters) designated as CEN22, CEN44 and CEN66 for second, fourth and sixth 

order respectively. 


The inviscid case uses a uniform Cartesian grid. The two compressible mixing layer test 
cases use a uniform Cartesian grid in the ^-direction and a mildly stretched Cartesian grid in the 
y-direction. In order to assess the true performance of the algorithm, no attempt is made to enhance 
the resolution using appropriate adaptive grid procedures. The code used for the Yee et aL( 19 ) 

study is employed for the present study. For non-periodic boundary conditions (BCs), the code 
reduces to lower order central base schemes near the boundary points. For the current study, we 
employ the same numerical BCs treatment in order to have a one-to-one comparison with the 
results obtained in Yee el al. (1999). Appropriate stable boundary difference operators developed 
by Strand (1994) should be used but are not yet implemented for the present study. The globa 
accuracy of the scheme related to intermediate BC treatment for the multi-stage Runge-Kutta 
method (Carpenter et al. ( 1995)) is not addressed here. Except for the vortex convection problem, 
all computations impose intermediate BC updates. Special treatment of the temporally (Carpenter 
et al ) and spatially dependent physical BCs related to the Runge-Kutta method ts not considered 
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for the two compressible mixing layer cases and such treatment is beyond the scope of this 
paper. Nonreflecting BCs or characteristic inflow and outflow boundary treatments are also not 
implemented. As indicated in the objective section, we explore the possible side benefits of 
the entropy splitting without considering the accompanying stable numerical boundary difference 
operator as a complete package for stability requirements. The three numerical examples were 
chosen to consist of periodic BCs, or computational domains whose boundaries are far enough 
away so as to not affect the mainstream flow activities. For the non-periodic cases, lower order 
non-characteristic boundary schemes are used. Evaluation of the performance of these schemes for 
the two compressible mixing layer test cases should take the above assumption into consideration. 


4.1. Isentropic Vortex Evolution 

The first test case is the evolution of a 2-D inviscid isentropic vortex in a free stream with 
periodic BCs in both spatial directions. The free-stream flow velocity, and v*,, pressure, 
Poo, and density, poo are (uoo,«<x>) = ( 1 , 0 ) andpoo = Poo = 1 An isentropic vortex with no 
perturbation in entropy ( SS = 0) is added to the free-stream flow field as initial conditions. The 
perturbation values are given by 


( 6u,Sv ) 



(~y,*)> 


(4.1.1) 


6T=- 


(7 ~ 1 ) 4 2 .i 

877T 2 



(4.1.2) 


where $ is the vortex strength and 7 = 1.4. Note that the vortex strength $ should not be confused 
with the /3 in section 2.1.2. Here T — Too = 1-0, (*,y) = ( z — *oo>y — y»o)> where *» 0 and 
y Vo are the initial coordinates of the center of the vortex, and r 2 — x 2 +y 2 . The entire flow field 
is required to be isentropic. Thus, for a perfect gas, = 1. 


From the relations, u = Uoo + Su, v = t>oo + Sv, T = Too + ST, and the isentropic relation, the 
resulting initial state for the conservative variables is given by 


p =T^ = (Too +ST)^ = 
pu =p(«oo + f> u ) = P 


!_ (7 — 


87 * 


l 

7-1 


8 1-* 3 

1 V 


(4.1.3) 

(4.1.4) 


pv =p(v oo + Sv) = p—e 1 x 


P =P y 


+ t;p( u2 + u *) 


7-1 2 


e 


(4.1.5) 

(4.1.6) 

(4.1.7) 
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The exact solution with given initial states is just a passive convection of the vortex with the 
free-stream velocity and thus provides a good measure of the accuracy of the schemes for smooth 
solutions of the nonlinear Euler equations. The vortex strength $ = 5 is fixed for all runs. Let k 
be the order of the central scheme; then the initial vortex covers a domain 0 < x < 10 + 0.125fc 
and -5 < y < 5 and its center is placed at (x Vo ,y Vo ) = (5,0). A uniform grid spacing of 
Aar = 0.125 and Ay = 10/(79 + k) is used. Although the actual grid size is 80 x 79, regardless 
of the order of the scheme, the grid size including ghost cells to accommodate the periodic BCs 
is (80 + k) x (79 + k). The reason for using an odd number of grid points in the y-direction is 
due to the compressible mixing layer structure of the code to accommodate fluctuations added to 
the inflow. The vortex is convected to the right by the mean flow velocity. Since there are no 
shock waves or steep gradient regions for this flow, the filter is used only to stabilize the nonlinear 
governing equations. For this reason, the filter coefficient k (2.4.2), if needed, should be kept very 
small. We use 0.001 < « < 0.07 for the computations. Due to the isentropic flow property, one 
can set p* (2.1.1 lc) to be a constant, p* = 1 is used for this test case. 

Density profiles at the centerline, y = 0, cutting through the center of the initial vortex are used 
for comparing the various schemes. Due to the time and spatial discretization numerical errors, 
the vortex, after long time integrations, can drift away from the centerline. The amount of drift 
depends on the scheme, grid size and the time step. If the computed vortex drifts away from the 
centerline but still preserves the vortex shape and strength, the centerline, y = 0, density profiles 
do not convey the full information and can be misleading. We complement the comparison with 
snap shots of density contours at different times up to 1 10 spatial periods. Here, one period is 
defined as the length of the periodic computational domain. The time required for one spatial 
period is t = 10. In all of the computations 6 = 0.1, where 6 is defined by (2.26a) of Yee et 
al. (1999). The limiter used is that given by Eqs. (2.250 of Yee et al. (1999). Here, 6 is the 
entropy satisfying parameter of Harten and Hyman ( 1983) for TVD schemes. As recommended by 
Carpenter et al. (1995), no intermediate BC update is imposed to improve the time accuracy of the 
multistage Runge-Kutta methods. The computations using intermediate BC updates do not have a 
drastic effect, but tend to diverge a little earlier than the case where there are no intermediate BC 
updates. 

We present results for sixth-order schemes. (Comparisons with second and fourth-order results 
are made in the following discussions.) Centerline distributions for CEN66, CEN66-ENT, ACM66, 
ACM66-ENT and ACM66S-ENT are shown in Figs. 4.1.1 - 4.1.7. Snap shots of selected density 
contours for these schemes are shown in Figs. 4. 1 .8 - 4. 1 . 16. 

The performance of the central schemes (with or without ACM) using entropy splitting and 
their un-split cousins is evaluated based primarily on vortex preservation capability after long time 
integrations of up to 130 periods (f = 1300). The discussion of numerical results is divided into 
the following: 

(a) Effect of the order of accuracy of the base scheme (At = 0.04): 

We ran some computations with second and fourth-order accurate base schemes for At 0.04 
(not shown), and compared them with the sixth-order results. CEN22 diverged after 3.5 periods, 
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CEN44 diverged after 6.8 periods, and CEN66 diverged after 5.17 periods. With no dissipation, 
the more accurate CEN66 computations had non-linear instability which caused it to diverge 
slightly earlier than the CEN44 case. Even so, the CEN66 density distribution was very accurate 
up to 5 periods, as seen from Fig. 4.1. Ic. It did show oscillations at period 5, indicating incipient 
instability at period 5.17. 

With the entropy splitting, the CEN22-ENT case diverged at 4.4 periods, the CEN44-ENT at 
13 periods, and the CEN66-ENT at 17 periods. The latter distribution is shown in Fig. 4.1.2c. 
Note that the entropy splitting allowed the calculation to proceed many periods further before it 
diverged, but, more significantly, the stability was improved as the order of accuracy was increased. 
This demonstrates the stabilizing effect of the entropy splitting. For a smaller A t = 0.01, the 
CEN66-ENT is stable up to 53 periods. See the later discussions. 

The addition of dissipation improved the performance markedly. Using the same At = 0.04, 
the ACM22 case with k = 0.07 diverged after 10 periods, the ACM44 with k = 0.06 became 
very distorted after 40 periods (with k, = 0.04 diverged after 32 periods), and the ACM66 with 
k = 0.06 remained stable for as long as we ran. (We stopped computing after 120 periods.) The 
density contours for this last case are shown in Fig. 4.1.12. Note that the vortex center starts to 
drift vertically and horizontally, and undergoes gradual smearing and distortion starting at period 
40 (see later discussion). The centerline distribution in Fig. 4.1.3c is therefore shown only up to 
period 30. The agreement with the exact solution is excellent. 

The addition of the entropy splitting yields further improvement. The ACM22-ENT with 
k = 0.07 remained stable, but the solution became very distorted at 30 periods. The ACM44-ENT 
with k = 0.04 similarly became severely distorted, and drifted vertically and horizontally beyond 
40 periods. The ACM66-ENT solution with k = 0.04 remained very good up to period 120, as 
shown in Fig. 4.1.13, although it started to drift to the right and upward even at period 30. This 
drift is evident from the centerline distributions in Fig. 4.1.4d. 

(b) Effect of ACM dissipation: 

We first will compare sixth-order computational results without dissipation (CEN66) with those 
due to added second-order upwind dissipation (ACM66). Results without entropy splitting for 
A t = 0.04,0.02, and 0.01 are shown in Figs. 4.1.1, 4.1.8 and 4.1.9 for CEN66, and Figs. 4.1.3 and 
4.1.12 for ACM66. The CEN66 computations diverged shortly after five periods for all three time 
steps. The ACM66 computations remained stable for as long as we ran (120 periods) for all three 
time steps. Due to the drift of the vortex center, the centerline distributions in Fig. 4.1.3 are shown 
only up to 30 periods. Note that up to that time, the larger time step, A t = 0.04 experienced less 
drift than the smaller time steps. The density contours for A t = 0.04 show smearing and distortion 
as we increase the duration of the time integration. 

Results with entropy splitting are shown in Figs. 4.1.2, 4.1.10 and 4.1. 1 1 for CEN66-ENT, and 
Figs. 4. 1 .4, 4. 1 .6, and 4. 1 . 13-4. 1.15 for ACM66-ENT. The entropy splitting for the non-dissipative 
computation (CEN66-ENT) allows a much longer time integration than that of CEN66 before it 
diverged, but even for the smallest time step, A t = 0.01, the solution for CEN66-ENT deteriorated 
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after 30 periods. Since ACM66-ENT allows a stable computation with less dissipation, the 
entropy splitting (« = .01 and At = 0.02,0.01) eliminated the distortion and smearing found in 
the ACM66 computations. The center of the vortex still drifted as we increased the time of the 
computation, but for the smallest time step. At = 0.01, there was only a very small drift to the 
right, even at 130 periods. Otherwise the vortex remained undistorted. Figures 4.1.5 and 4.1.6 
show the centerline density distribution for 10 — 130 periods with a 10 period increment. One can 
see the drifting effect as a function of the At. Aside from the drifting, the vortex is still quite 
accurate. This is evident from the density contours Figs. 4.1.14 and 4. 1.15. 


We would like to point out that the vertical drifting of the vortex away from the centerline 
y - o and horizontal drifting (or rather shifting) are quite common for all schemes beyond 30 
periods. Depending on the scheme, the amount of numerical dissipation and the time step, drifting 
can occur as early as 5 periods. We believe that the vertical drifting is due largely to the spatial 
numerical dissipation of the scheme. This is evident from the k refinement study. See Figs. 4.1.13 
and 4.1.15. The horizontal drifting is due largely to the phase error of the time integrator. This 
is evident from the time step refinement study on ACM66-ENT using « = 0.01. See Figs. 4.1.5, 
4.1.6, 4.1.14 and 4.1. 15. 

(c) Effect of the adjustable ACM parameter k : 

Although we experimented with various values of k (0.01 < k < 0.07) to find the optimum 
value, we show some results only for the case giving the best solution, namely ACM66-ENT. 
For the larger time step. At = 0.04, Fig. 4.1.4 shows a negligible effect up to 30 periods when 
increasing the value of k from 0.01 (Fig. 4.1.4c) to 0.04 (Fig. 4.1.4d). The computation k = 0.04 
and At = 0.04 remained stable, although with some distortion and smearing, as shown by the 
density contours in Fig. 4.1.13. Actually, this case gives slightly better results than the one using 
k = 0.01 and the same time step size (not shown). The computation is not stable for k < 0.01 and 
At > 0.01 for the ACM66-ENT scheme. 

(d) Effect of entropy splitting: 

The cases discussed above clearly show the advantages of using entropy splitting. For the cases 
without dissipation, CEN66-ENT, the splitting allowed the computation to proceed for a much 
longer time before it became unstable. The real advantage came when used in conjunction with 
the upwind filter, ACM66-ENT, where excellent solutions with just a very small amount of drift 
were obtained after long time integrations of 130 periods. 

(e) Effect of the time step: 

In general, decreasing the time step At gave a better solution. For the base scheme CEN66, 
the improvement in decreasing the time step from 0.04 to 0.01 was negligible. However, with 
entropy splitting, decreasing the time step had a significant effect, as seen from the results for 
CEN66-ENT in Figs. 4.1.2, 4.1.11 and 4.1.10. The improvement with decreasing time step was 
not as marked when the ACM filter was applied. ( Fig. 4.1.3 showing the density distributions for 
ACM66 actually appears to show an improvement with increasing time step. This is misleading, 
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since the slight upward drift of the vortex center produced a shift in the centerline distributions for 
the smaller time steps. The density contours illustrate similar vortex preservation capability (not 
shown)) The improvement in decreasing the time step from 0.02 to 0.01 is clearly shown in the 
density distribution and density contours for ACM66-ENT of Figs. 4.1.5, 4.1.6, 4.1.14 and 4. 1.15. 
The significant downward drift for At = 0.02 has been totally eliminated for At = 0.01. The 
drift to the right has also been virtually eliminated. A quantitative evaluation of the solution for 
At = 0.01 can be obtained from the centerline distribution for ACM66-ENT in Fig. 4.1.6. Even 
at 130 periods, the profile is undistorted, with a shift due to the small drift to the right. 

(f) Effect of symmetric vs. upwind ACM: 

We confine the comparison to the best solution, namely ACM66-ENT for At = 0.01. The 
centerline distributions up to period 30 for the symmetric ACM case, ACM66S-ENT in Fig. 4.1.7 
are as good as for the upwind case, ACM66-ENT in Fig. 4. 1.4a. The symmetric TVD dissipation of 
Yee (1985), limiter 2.7b was used. Note that the value of k in the former case has been decreased 
to 0.005. A more meaningful comparison can be made by examining the density contours up to 
1 10 periods for the two cases in Figs. 4.1.16 and 4.1.15. Note that the symmetric ACM solution 
undergoes some distortion and upward drift as the period is increased. Actually ACM66S-ENT 
with k = 0.005 and At = 0.01 produce better results than ACM66-ENT with k = 0.04 and 
At = 0.01 or k = 0.01 and At = 0.04 The drifting behavior also occurs with ACM66-ENT if 
k > 0.04 or for k = 0.01 and At = 0.04 If we had used a larger value of k = 0.01 for the 
symmetric case, the solution would have become very inaccurate. On the other hand, a smaller 
value of k = .001 produced an unstable solution due to insufficient dissipation. It appears that 
the use of symmetric ACM filter in conjunction with entropy splitting is also computationally 
attractive. In all of the computations using ACM, the solution is quite sensitive to the value of k 
and the time step size, although, the upwind ACM appears to be a bit less sensitive. 


(g) Effect of the splitting parameter /3: 

All of the calculations shown for the entropy splitting have been for a value of /3 = 1 
(a = —1.8). This produces an equal amount of conservative and non-conservative splitting. 
We have also run some cases for/3 = 0.5 (a = -1.6). This gives a splitting that is one-third 
conservative and two-thirds non-conservative. The results are slightly worse than for the (3 = 1 
case. Increasing the conservative proportion beyond 80% will defeat the purpose of using the 
splitting for this particular example since the gain in stability is diminished by the expense of the 
added CPU computation required by the splitting. 

In summary, the use of entropy splitting in conjunction with an upwind TVD ACM filter has 
preserved a horizontally convecting vortex with great accuracy after long time integration of 130 
periods. The splitting helps minimize the use of numerical dissipation. To the authors knowledge, 
highly accurate finite discretization computations previously reported in the literature were only 
carried out up to 10 periods of integration. 
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4.2. Vortex Pairing in a Time- Developing Mixing Layer 

This test case involved vortex growth and pairing in a temporal mixing layer at a convective 
Mach number equal to 0.8. At this Mach number there are shock waves (shocklets) that form 
around the vortices and the problem is to compute accurately the vortex evolution while avoiding 
oscillations around the shocks. Previous calculations of the problem can be found in Sandham 
and Reynolds (1989), Lumpp (1996), Fu and Ma (1997), Sandham and Yee (1998) and Yee et al. 
(1999). Figure 4.2.1 shows a schematic of the physical problem. Here we set up a base flow as in 
Sandham and Yee (1989) 


u = 0.5 tanh(2y), (4.2.1) 

with velocities normalized by the velocity jump v,\ — u 2 across the shear layer and distances 
normalized by the vorticity thickness. 


■ - Ul ~ u » 

( du/dy)mam 


(4.2.2) 


Subscripts 1 and 2 refer to the upper (y > 0) and lower ( y < 0) streams of fluid respectively. The 
normalized temperature and hence local sound speed squared is determined from an assumption of 
constant stagnation enthalpy 


c 1 = c\ + - »'). (4.2.3) 

Equal pressure through the mixing layer is assumed. Therefore, for this configuration of u 2 = -«i 
both fluid streams have the same density and temperature for y —* ±oo . The Reynolds number 
defined by the velocity jump, vorticity thickness and kinematic viscosity at the free-stream 
temperature is set equal to 1000. The Prandtl number is set to 0.72, the ratio of specific heats 
is taken as 7 = 1.4, and Sutherland’s law with reference temperature Tr = 300 o iif is used for 
the viscosity variation with temperature. The reference sound speed squared, c\, is taken as the 
average of c 2 over the two free streams. 

Disturbances are added to the velocity components in the form of simple waves. For the normal 
component of velocity we have the perturbation 

2 

v' = ^ an cos(2nkx/ L a + 4>k)cxp(-y 2 /b), (4-2.4) 

A=1 


where L a = 30 is the box length in the ^-direction and b = 10 is the y- modulation. In our test 
case we simulate pairing in the center of the computational box, by choosing the initially most 
unstable wave k = 2 to have amplitude a 2 = 0.05 and phase fa = —n/2, and the subharmonic 
wave k = 1 with a x = 0.01 and fa = -tt/ 2. The u-velocity perturbations are found by assuming 
that the total perturbation is divergence free. Even though these fluctuations correspond only 
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approximately to eigenfunctions of the linear stability problem for a compressible mixing layer, 
they serve the purpose of initiating the instability of the mixing layer and have the advantage as a 
test case in that they can be easily coded. 

Numerically the grid is equally spaced and periodic in the x — direction and stretched in the 
y- direction, using the mapping 


L y sinh( 6 y J 7 ) 
2 sinh( 6 y ) ’ 


(4.2.5) 


where we take the box size in the y-direction Ly = 100, and the stretching factor b y = 3.4. The 
mapped coordinate tj is equally spaced and runs from —1 to +1. The boundaries at ±2/ y /2 are 
taken to be slip walls. For example, at the lower boundary 


Pi = Pi-, (4.2.6a) 

(pu) i = {pu) 2 , (4.2.6b) 

(pv) i = 0, (4.2.6c) 

(e)i = [4(e) j - (e),]/3, (4.2.6d) 


where subscripts here refer to the grid point and e is the total energy. 


We compute this test case on a 101 x 101 grid. A gird refinement study was performed in Yee 
et al. (1999) Figure 4.2.2, a reference solution taken from Yee et al. (1999) using ACM44, shows 
a snapshot of the temperature contours at t = 40, 80, 120 and 160 using ACM44, illustrating the 
roll-up of the primary vortices followed by vortex merging. Shock waves and shears form around 
the vortices with a peak Mach number ahead of the vortex of approximately 1.55 at t = 120. The 
grid is 201 x 201 . 

For this vortex pairing in a time-developing mixing layer, we study only the effect of the choice 
of the arbitrary splitting parameter /3 (i.e., the proportion of conservative and non-conservative 
parts of the splitting) in obtaining the same shock location as the un-split approach with scheme 
ACM 66 -ENT using Af = 0.1. In all of the computations for the vortex pairing case, limiter 
(2.25h), and 6 = 0.25 (2.26a) of Yee et al. (1999) are used and k = 0.7 (2.4.2) is used for the 
nonlinear fields for the ACM methods. Intermediate BC update is imposed in order to have a 
one-to-one comparison with the Yee et al. (1999) results. 

We consider a — — 100 , - 10 , - 5 , —3, —2, -1.8 (/3 = 1), —1.6, 0.1, 1 (/3 = — 6 ), 
2, 5, 10, 100 with /3 = (a + i)/(l - 7 ). The scheme diverges for a = 0.1. This corre- 
sponds to 136.36% of the conservative proportion and -36.36% the non-conservative proportion. 

By monitoring the left A-shock location, studies indicate that for a > |5|, the same shock 
location and shock strength of the A-shock are obtained as in the un-split approach. With the 
exception of a small increase in spurious noise in the vicinity of the shock (not shown), it is 
surprising to see that a slightly over 100 % conservative proportion (a > 0 ) and the corresponding 
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negative non-conservative proportion would give the correct shock strength and location. Away 
from the A-shock area, the solution is less sensitive to the choice of a > |5|. 

For the physical choice of a < - 7 , we obtain the opposite effect as compared to that of a > 0, 
in terms of spurious noise. As a ^ 1 - 27 , the entropy splitting has a spurious noise reduction 
capability when compared with the un-split approach. For example, the shock strength and location 
that are a bit away from the A-shock are almost the same for a = -3 (/3 = 4) as for a = -5 
0 _ 9 ) except that a = -3 has a bigger smoothing effect on the spurious noise generated 
by the scheme (especially when a more compressive flux limiter is employed). In addition, for 
_5 < a < — y • bigger negative a in that range results in more shift of the A-shock location away 
from the un-split approach location. For example, for a = —3, there is a shift of approximately 
1_1 1/2 grid points. The shock strength reduction at the A-shock location is very small. A 
reduction in At might improve the accuracy of the shift in the shock location, based on the vortex 
convection case. Here, we only compare the results with the un-split approach using the same time 
step and BCs as reported in Yee et al. (1999). Figure 4.2.3 illustrates the normalized temperature 
at t = 160 using a = -5 and a = -3 on a 101 x 101 grid. Except for a slight noise reduction in 
the vicinity of the shock, the ct = — 5 solution is almost identical to the un-split computation. To 
illustrate the noise reduction capability of the splitting. Fig. 4.2.4 shows a comparison of the split 
and un-split forms with the ACM filter turned off for the linearly degenerate fields (u,u and v,v 
characteristic fields) using a = —3 (J3 = 4). One can see the noise reduction effect of the entropy 
splitting on the scheme which, at the same time, maintains the accuracy of the shock and shears 
away from the A-shock location as in the un-split approach. Since the entropy splitting requires the 
same amount of filter as the un-split approach for this type of rapidly developing shock-turbulence 
interactions, its stabilizing effect is only on the spurious noise reduction, and the benefit is not 
as pronounced as for the smooth flow case. However, for turbulent flows involving long time 
integrations that contain weak or no shock waves, the entropy splitting could help minimize the 
use of numerical dissipation due to its unique nonlinear stability property. A separate investigation 
is in progress. 

4.3. Shock Wave Impingement on a Spatially-Evolving Mixing Layer 

The third test case has been developed to test the behavior of the schemes for shock waves 
interacting with shear layers where the vortices arising from shear layer instability are forced to 
pass through a shock wave. Figure 4.3.1 shows the schematic of the physical problem. An oblique 
shock is made to impact on a spatially-developing mixing layer at an initial convective Mach 
number of 0.6. The shear layer vortices pass through the shock system and later through another 
shock, imposed by reflection from a (slip) wall at the lower boundary. The problem has been 
arranged with the Mach number at the outflow boundary everywhere supersonic so that no explicit 
outflow boundary conditions are required. This allows us to focus on properties of the numerical 
schemes rather than on the boundary treatment. 

Figure 4.3.2, a reference solution taken from Yee et al.( 1999) using ACM44, shows the nature of 
the flow on a 641 x 161 grid illustrating the pressure, density and temperature fields using the 
ACM44 method and the same limiter as the pairing case with k = 0.35 for nonlinear characteristic 
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fields and k = 0.175 for linear characteristic fields. The time step is At = 0.12. The entropy 
satisfying parameter 6 of Harten and Hyman (1983) is set to 0.25. An oblique shock originates 
from the top left hand comer and this impacts on the shear layer at around x = 90. The shear layer 
is deflected by the interaction. Afterwards we have a shock wave below the shear layer and an 
expansion fan above it. The shock wave reflects from the lower solid wall and passes back through 
the shear layer. The lower wall uses a slip condition so no viscous boundary layer forms and we 
focus on the shock-wave interaction with the unstable shear layer. 

The inflow is specified again with a hyperbolic tangent profile, this time as 


u = 2.5 + 0.5 tanh(2y), (4.3.1) 

giving a mixing layer with upper velocity u 2 = 3, lower velocity u 2 = 2, and hence a velocity 
ratio of 1.5. Equal pressures and stagnation enthalpies are assumed for the two streams, with 
convective Mach number, defined by 


M c = 


«1 ~ 

Cl + C 2 ’ 


(4.3.2) 


where c x and c 2 are the free stream sound. speeds, equal to 0.6. The reference density is taken as 
the average of the two free streams and a reference pressure as ( pi + /> 2 )(u i — u 2 ) 2 /2. This allows 
one to compute the inflow parameters as given in the first two columns of Table 4.3.1. Inflow 
sound speed squared is found from the relation for constant stagnation enthalpy (4.2.3). The 9 in 
Table 4.3.1 is the flow inclination angle with respect to the ^-direction. 

The upper boundary condition given in column 3 of Table 4.3.1 is taken from the flow properties 
behind an oblique shock with angle /3 = 12°, The table also gives the properties behind the 
expansion fan (column 4) and after the oblique shock on the lower stream of fluid (column 5) 
computed by standard gasdynamics methods with /3 = 38.118°. In practice, the conditions in 
regions 4 and 5 do not correspond exactly to the simulations due to the finite thickness of the 
shear layer. The Mach number of the lower stream after this shock is approximately M 5 = 1.6335 
and remains supersonic through all the successive shocks and expansion fans up to the outflow 
boundary. The resulting shock waves are not strong, but tests showed that they could not be 
computed without using shock-capturing techniques. The lower boundary was specified with the 
same slip condition used for the pairing case (Equation (4.2.6)). 

The Prandtl number and ratio of specific heats were taken to be the same as for the vortex pairing 
test case. The Reynolds number was chosen to be 500. 

Fluctuations are added to the inflow as 


v' = cos(2nkt/T + <f>k)exp(-y 2 /b), 

h=l 


(4.3.3) 
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with period T = A/u c , wavelength A = 30, convective velocity u c = 2.68 (defined by 
u c = (uic 2 + u 2 c 1 )/(c 1 + c 2 )) and 6 = 10. For* = 1 we take a 2 = 0.05 and <j> = 0, and for fc = 2 
we take a 2 = 0.05 and <f> = jt/2. No perturbations are added to the u-component of velocity. 

The grid is taken to be uniform in * and stretched in y according to equation (4.2.5) with b y = 1. 
This stretching is much milder than for the pairing problem, as we have to resolve the shear layer 
even when it deflects away from y = 0. The box lengths were taken to be L, = 200 and L y = 40. 

The reference solution indicates that vortex cores are located by low pressure regions and the 
stagnation zones between vortices by high pressure regions. The shock waves are seen to be 
deformed by the passage of the vortices. Another interesting observation is the way the core of 
the vortex at * = 148 has been split into two by its passage through the reflected shock wave. 
In spite of the relatively high amplitude chosen for the subharmonic inflow perturbation, there is 
no pairing of vortices within the computational box. We do, however, begin to see eddy shock 
waves around the vortices near the end of the computational box where the local convective Mach 
number has increased to around 0.66. The oscillations seen near the upper boundary for z > 120 
occur where the small Mach waves from the initial perturbations arrive at the upper boundary. The 
use of characteristic boundary conditions should remove this problem. Practically, the amplitude 
of oscillations is not sufficient to cause numerical instability or affect the remainder of the flow. 

For this shock wave impingement on a spatially-evolving mixing layer, again, we study 
only the effect of the choice of the arbitrary splitting parameter 0 in obtaining the same 
shock location as the un-split approach. The study is limited to ACM 66 -ENT and ACM 66 
using the fifth limiter of equation (2.25h) of Yee et al. (1999), and the same k value and 
time step size (At = 0.12) as the reference solution. Intermediate BC updates are imposed 
in order to have a one-to-one comparison with the Yee et al. (1999) result. We consider 
a = - 100 , - 10 , -5, - 3 , - 2 , - 1.8 ((3 = 1), - 1 . 6 , 0.1, 1 ((3 = - 6 ), 2 , 3, 4, 5, 10 , 20, 100. 
Studies indicate that for a = -2, -1.8, -1.6, 0.1, 1, 2, 3, 4, 5 the solution diverges. The rest of 
the a. values for |ct| > 10 produce almost identical results as the un-split case. This example poses 
a more stringent requirement on the a range than the vortex pairing case. Figure 4.3.3 compares 
the un-split pressure contours with the split case for a = rfclO at t = 120 on a 321 x 81 grid with 
At = 0.12. The value of a = +10 (J3 = -28.5) corresponds to a 103.6% conservative proportion 
and a small negative non-conservative proportion. The a = — 10 (J3 = 21.5) corresponds to a 
less than 100% conservative proportion. Note that the a = -10 solution produces spurious noise 
reduction, while the ol — +10 solution actually induces more spurious noise than the un-split 
approach. This opposite effect of the spurious noise phenomena for a > 0 and a < -7 is shared 
with the vortex pairing example. Again, for turbulent flows involving long time integrations that 
contain weak or no shock waves, we believe the entropy splitting could help minimize the use of 
numerical dissipation due to its unique nonlinear stability property. 

4.4. Computational Costs 

For the compressible mixing layer computations using the fourth-order Runge-Kutta method, 
the central base schemes with the ACM filter are only around 25% more expensive than the same 
base schemes without ACM filter. This has been achieved by only requiring one application of 
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the ACM filter per full time step for the convection terms. For LMM time discretizations, the 
central base scheme with the ACM filter is only 10% more expensive than standard second-order 
TVD schemes. The entropy splitting is approximately 20 % more expensive than the un-split 
conservative form for the 2-D mixing layer computations in conjunction with the fourth-order 
Runge-Kutta method. The extra CPU time is mainly due to the fact that, for each direction, four 
entropy splittings are required. If two to three time level LMM type of time discretizations are 
used, less CPU time can be realized. 


Summary 


Our study shows that the entropy splitting can be formally extended to a thermally perfect gas, 
with the internal energy being an arbitrary function of temperature. For non-equilibrium flows 
consist of a mixture of different species, each obeying a thermally perfect gas law, extension of the 
splitting is problematic. While we were able to prove the symmetry and homogeneity properties, 
the degree of homogeneity can only be obtained by solving a system of nonlinear equations. 
In addition, to obtain the Jacobian of the transformation required inverting a non-sparse linear 
system. It would therefore be difficult to establish the positive definite condition. Consequently, 
the extension of the method to non-equilibrium flows is not practically feasible. If the homogeneity 
condition is not required, then one can use symmetry variables based on the physical entropy, as 
was shown by Chalot et al. (1990). In this case, the resulting PDEs are in pure non-conservative 
form and entropy splitting is no longer possible. For magnetohydrodynamics, the magnetic field 
has to be added as a "conservative" variable. But the square of the magnetic field is one of the 
terms in the definition of the total energy. Thus, from dimensional arguments it is clear that one 
cannot obtain the homogeneity condition. A similar situation exists for the artificial compressibility 
method of treating incompressible flow. 

Using the same high order central schemes, numerical experiments with a 2-D vortex convection 
Euler computation consisting of periodic BCs indicate that entropy splitting is more stable than the 
un-split (purely conservative) approach. With an appropriate time step, numerical dissipation is 
not required for up to 30 spatial periods with nearly perfect vortex preservation as opposed to only 
5 periods for its un-split cousin. For even longer time integration, although numerical dissipation 
is needed to stabilize the schemes, the amount required is much less than for its un-split cousin. A 
nearly perfect vortex preservation of up to 120 periods was achieved. 

For the mixing layer study, in order to obtain the same A-shock strength and shock location 
as the un-split approach using the same scheme, the range of the arbitrary splitting parameter 0 
has to be confined to the use of at least 90% of the conservative proportion of the flux derivative. 
For problems without A-shocks, a wider range of (3 can be used. Only a slight advantage of 
the entropy splitting over the un-split approach was observed for this type of flow physics. The 
advantage is in terms of noise reduction and improved nonlinear stability. There is no reduction in 
the use of the ACM filter. This might largely be due to the rapidly developing flow and the high 
percentage of conservative proportion required. Unlike the vortex convection with periodic BCs 
in all spatial directions, these two more complicated cases consist of rapidly developing Hows. 
Not all of the physical BCs are periodic. In addition, the BCs consist of spatially or temporally 
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sinusoidal disturbances. Thus, the performance of the schemes for the mixing layercases is partially 
clouded by the spatially and temporally dependent physical BCs that required special treatment 
in conjunction with the Runge-Kutta time integrator, and also by the fact that we did not impose 
the more stable and appropriate boundary difference operator. The use of the symmetric form of 
Harten (1983a) for the viscous term might be a source of improvement. In addition, a time step 
reduction will also benefit the use of entropy splitting as indicated in the vortex convection study. 
Without additional study, the benefit of using the entropy splitting is inconclusive for compressible 
turbulence mixing applications. Since the entropy splitting requires the same amount of filter as 
the un-split approach for rapidly developing shock-turbulence interactions, its stabilizing effect 
is not as pronounced as for the smooth flow case. However, for turbulent flows involving long 
time integrations that contain weak or no shock waves, the entropy splitting could help minimize 
the use of numerical dissipation due to its unique nonlinear stability property. More rigorous 
implementation of the BCs and extensive study are needed. These are the subjects of future 
research. 

Overall, the three numerical examples indicate a positive side benefit of the entropy splitting. 
The splitting can stabilize spurious noise generated by the non-dissipative or low dissipative spatial 
discretizations which are a major cause of nonlinear instability. Modem high-resolution numerical 
dissipation has been the major player in improving nonlinear instabilities for short or moderate 
time integrations (unsteady). Most often, added numerical dissipation is necessary for longer 
time integration at the expense of excess smearing of the flow physics without resorting to finer 
grids and extremely small time steps. The use of the entropy splitting form of the flux derivative 
fn conjunction with high-resolution filters can minimize the use of numerical dissipation. We 
believe that the use of the entropy splitting is not limited to spatial central schemes (compact or 
non-compact), but to spectral and spectral-like spatial schemes as well. This and the blending of 
ACM filters with other filters on the possible suppression of spurious high frequency oscillations 
as discussed in section 2.6 are the subjects of our near future research. 
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APPENDIX A 

Entropy Splitting for a Thermally Perfect Gas in 2-D Cartesian Coordinates 


In this appendix we specialize the equations derived in Sections 3.1 and 3.2 to 2-D, stationary 
Cartesian coordinates. The variables in (2.1.1) are defined as 


" p m 


m 


n 

m 

n 

, F = 

m 2 / p -1- p 

mv 

, <7 = 

nu 

n 2 Jp + p 

. e . 


_(e +p)m/p. 


_(e +p)n/p. 


where m — pu, n — pv. u and v are the Cartesian velocity components, and the other variables 


41 


are as defined in Section 3. 1 . The temperature T(U) is obtained by solving implicitly the equation 

(A.2) 




P 2 p 2 

The transformed variable W takes the form 

W = [ w x w 2 Wi Wtf =-[e-2Z-p(l+/3) —pu -pv p j , 

P 


(A.3) 


where { is given by (3.1.19b). Note that wi corresponds to the component uf, w 3 andu>j correspond 
to the components w, and w 4 corresponds to the component w in Section 3. 1 . The upper triangular 
part of Uw is now written as 


Uw = j 


ap apu apv 

apu 2 -p apuv 

.2 


ae + bp 

u[ae + (6 - l)p] 
apv" — p v[ae + (b — l)p] 


(A.4) 


where q 2 = u 2 + t? 2 and a and b are given by (3.1.21b) and (3.1.21c). The upper triangular parts 
of Fw and Gw can be written as 


Fw = j 


apu apu 2 — p apuv u [ae + (b — l)p] 

u(apu 2 — 3p) v(apu 2 - p) u 2 {ae + (b - 2)p] - £(e + p) 
u(apv 2 — p) uv[ae + (b — 2 )p] 

/44 


(A.5) 


where 


and 


ae 2 _ x e 2, . P 2 


Gw = j 


hi. = «[— + P{2(6 - 1 )- - q 2 } + - WI +0)~ 2 }]> 
P P P 


apv apuv apv 2 — p w[ae + (& — l)p] 

v(apu 2 — p) u(apv 2 — p) uv[ae + (6 — 2)p] 

v(apv 2 — 3 p) v 2 [ae + (b — 2)p] — jj(e + p) 

<744 


(A.6) 


(A.7) 


where 


,ae 2 . p 2 


<744 = V[— +p{2(6 - 1)- - q 2 } + ^{Kl+0) ~ 2}]. 
P P P 


The matrix A a = Fu can be written as 


A x = 


0 10 0 
K\ — u 2 (2 — k)u —kv k 

—uv v u 0 

(Ki — H)u H — KU 2 —KUV (1 + k)u 


(A.8) 


(A.9) 
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where K\ — -f- H ~ h J9 2 is the total enthalpy per unit mass, h — c -j~ T is the specific 

enthalpy, and \ and K are gi ven by (3.2.2). The three distinct eigenvalues of A a are 

A, 1 = u, A , 2 =u + c, and A, 1 = u - c, (A. 10) 


where c is given by (3.2.5). The right eigenvector matrix R„ can be written as 


R, = 


r l o 

u o 

V c 


1 

u + c 

V 


1 

u — c 

V 


lk 2 vc H+uc H — uc J 


(A.ll) 


where K 2 = - x/k, and the first two columns correspond to A, 1 . The column vector 

R^AU is given by 


R a ~ l AU = 


A p - Ap/c 2 

pAv/c 

\{Ap/c 2 + pAu/c) 
[|(Ap/c 2 — pAu/c) J 

The Roe-averaged states are defined in Section 3.2. 

The matrix A y = Gu can be written as 



0 

0 

1 0 - 


— uv 

V 

u 0 

II 

Kt - v 2 

— KU 

(2 — k)v k 


.(Kt-H)v 

-/cur 

H - KV 2 (1 + k)u- 

The three distinct eigenvalues of A y are 



Ay 1 = V , Ay 2 = 

V + c, 

and Ay 3 = v — c. 


The right eigenvector matrix R y can be written as 


Ry = 


r l 

u 

V 


0 

c 

0 


1 

u 

V + c 


1 

u 

v — c 


L k 2 uc H + vc H — vc J 


(A. 12) 


(A.13) 


(A. 14) 


(A.15) 


The column vector R y 1 AU is given by 


R y ~ 1 AU = 


A p - Ap/c 2 

pAu/c 

j(A p/c 2 + pAv/c) 

IWAp/c' -pA./S) 


(A.16) 
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If we set e = 1/(7 - 1), and make further simplifications, (A. 3 - A. 16) reduce to the perfect 
gas case. 


APPENDIX B 

Entropy Splitting for a Thermally Perfect Gas in 3D Cartesian Coordinates 

In this appendix we specialize the equations derived in Sections 3.1 and 3.2 to 3D, stationary 
Cartesian coordinates. The conservation law in this case can be written in the form 

U t + E, + F y + G t = 0, (B.l) 

where U t = ^,E. = ff ,F y = and G z = ^. The dependent variable U is the vector of 
conservative variables, while E, F and G are the flux vector components. The variables in (B.l) 
are defined as 


U = 


' p " 
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n 2 /p +p 
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.(e + p)l/p. 


_(e +p)m/p. 


.(e +p)n/p. 


(B.2) 


where l — pu, tti = pv . n — pw. u, v and w are the Cartesian velocity components, and the other 
variables are as defined in Section 3.1. The temperature T(U) is obtained by solving implicitly the 
equation 

e 1(/ 2 +m 2 +n 2 ) 

) ~ " ~ o „2 


(B.3) 


The transformed variable W takes the form 


W — [wi w-i utj ty 4 


ty 5 ] T = -fe - 2e — p(l + fi) -pu -pv -pw 
P 


(B.4) 


apu 


U w = -. 


apuv 
apv 2 — p 


apuw 
apvw 
apw 2 — p 


where ^ is given by (3.1.19b). Note that w 2 , and correspond to the components of w in 
Section 3.1 and w$ corresponds to w in Section 3.1. The upper triangular part of U\v is now 
written as 

' ap apu apv apw ae + bp 

u[ae + (b - l)p] 

v[ae + (b - l)p] | , (B.5) 

w[ae + (b - l)p] 

^+p(7 £ -5 1 )+^(1+«J 

where q 2 = u 2 + v 2 + w 2 and a and b are given by (3.1.21b) and (3.1.21c). The upper triangular 
parts of Ew Fw and Gw can be written as 

■ apu apu 2 — p apuv apuw u[ae + (b - l)p] 

j u(apu 2 — 3p) t>(apti 2 — p) iu(apti 2 - p) e 2 s 

Ew = ~ u(apv 2 — p) apuvw uv[ae -f (b — 2)p] | , (B.6) 

^ u(apw 2 — p) uw[ae + (b — 2)p] 
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where 


and 


e 2 s = u 2 [ae + (6 - 2)p] (e + p), 

/> 


e „ = u[— + p{2(b - 1)- - q 2 } + ^{6(1 + 0) - 2}]; 
P P P 


(B.7) 

(B.8) 


•Fw = ^ 


where 


and 


and 


apv apuv apv 2 — p apvw v[ae + (b- l)p] ] 

v(apu 2 — p) u(apv 2 — p) apuvw uv(ae + (6 — 2)p] 

v(apv 2 — 3p) w(apv 2 — p) fzs 

v(apw 2 — p) vw[ae + (b — 2)p] 

/s 5 


/j5 = u J [ae + (6 - 2)p] - -(e +p), 


/s5 = »[— - + p{2(6 - 1) — g 2 } + — (K 1 + 0) - 2 }]; 

P P P 


, (B.9) 


(B.10) 

(B.ll) 


Gw = j 


apw apuw apvw apw 2 — p w[ae + (b~ l)p] 

w(apu 2 — p) apuvw u(apw 2 — p) uw[ae + (6 — 2)p] 

w[apv 2 — p ) v(apw 2 — p) vw[ae 4- {b — 2)p] 

.2 


where 


and 


w(apw 2 — 3 p) 


045 = u> 2 [ae + (6 - 2)p] - — (e + p)> 


055 = ip[— - + p{2(& — 1) g 2 } H {&(1 + 0) — 2}]. 

p P P 


045 

055 


(B-12) 

(B.13) 

(B.14) 


The matrix A, = Eu can be written as 


A, = 


0 10 0 
K\ — u 2 (2 — k)u —kv —kw 
— uv v u 0 

—UW W 0 u 


0 

K 

0 

0 


[.(#!- ff)u H — KU 2 —KUV —KUW (1 + k)uJ 


(B.15) 


where K x = \nq 2 + *, H = h + |g 3 is the total enthalpy per unit mass, h = e + T is the specific 
enthalpy, and x and « are given by (3.2.2). The three distinct eigenvalues of A m are 


A, 1 = u, A a 2 - u + c, and A, J = u - c, 


(B.16) 
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where c is given by (3.2.5). The right eigenvector matrix R m can be written as 


- 1 

0 

0 

1 

1 - 

u 

0 

0 

u -j- c 

u — c 

V 

c 

0 

V 

V 

w 

0 

c 

w 

w 

-K 2 

VC 

wc 

l l - hue 

H - uc. 


where K 2 = - x/«, and the first three columns correspond to A, 1 . The 

R a ~ 1 AU is given by 


R,~ l AU = 


A p — A p/c 2 

pAv/c 

pAw/c 


|(Ap/c 2 + pAu/c) 

L|(Ap/c 2 — pAu/c) J 


The Roe-averaged states are defined in Section 3.2. 


The matrix A y = Fu can be written as 




0 

0 

l 

0 

—uv 

t; 

u 

0 

Ki - v 2 

— KU 

(2 - k)u 

— KW 

— wv 

0 

w 

V 

H 

1 

3 

— KUV 

H - KV 2 

— KWV 


0 

0 


K 

0 

(1 -I- k)v. 


The three distinct eigenvalues of A y are 


Ay 1 = v, Ay 2 = v + c, and A y * = v - c. 


The right eigenvector matrix R y can be written as 



■ 1 0 

0 1 1 • 


u c 

0 u u 

Ry — 

V 0 

0 v + c v — c 


w 0 

C ID W 


. K 2 uc 

wc H + vc H — vc. 

The column vector R y 1 AU is given by 


Ry 

II 

t=> 

<J 

A p — A p/c 2 

pAu/c 

pAw/c 

|(A p/c 2 + pAv/c) 

■ |(Ap/c 2 - pAv/c). 


(B.17) 
column vector 


(B.18) 


(B.19) 


(B.20) 


(B.21) 


(B.22) 
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The matrix A z = Gu can be written as 

0 0 


A z = 


—uw w 

—vw 0 

K\ - w i —ku 

(K 1 — H)w —KUXJD 


0 

0 

w 


1 

u 

V 

(2 — k)u> 


— KV 

— KVW H — KW 2 


0 
0 
0 
K 

(1 + k)w J 


The three distinct eigenvalues of A z are 


A* 1 = w 

II 

N 

w + c, and A x 3 = 

w — 

The right eigenvector matrix R z 

can be written as 



■ 1 0 

0 1 1 

- 


u c 

0 u u 


Rz = 

V 0 

c V V 



w 0 

0 w + c w — 

c 


. K 2 uc 

vc H + wc H — wc J 

The column vector R z 1 A U is given by 





A p — A p/c 2 




pAu/c 


Rz 

U = 

pAv/c 

. 



|(A p/c 2 + pAw/c) 




-|(Ap/c J -pAw/c). 



(B.23) 


(B.24) 


(B.25) 


(B.26) 


APPENDIX C 

Caloric Equation for an Ideal Diatomic Gas 

In this appendix we present the internal energy functions for an ideal diatomic gas. This gas 
is one for which rotation is fully excited, the vibrating molecule is approximated as a harmonic 
oscillator, and electronic contribution is neglected. The harmonic oscillator model implies that 
the temperature is small compared to the dissociation temperature. Under these assumptions, the 
internal energy is given by 

e(T) = 2 T + FTI' (C1) 

where 6 is the characteristic temperature for vibration, and 

r = BIT. (0.2) 

The quantity k = dt/dT is needed to solve Eq. (3.1.8) iteratively, and to evaluate the functions a, 
b, and \ in (3.1.2lb), (3.1.21c), and (3.2.2). It is given by 


rV 


* r) " 2 + 


(C.3) 
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The function f(T) defined in (3. 1 .6) is given by 


/(f) =?-•/»(.' -i)«p- 


re 


r 


e r - 1 


(C.4) 
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Property 

(1) 

(2) 

(3) 

(4) 

(5) 

u-velocity 

3.0000 

2.0000 

2.9709 

2.9792 

1.9001 

v - velocity 

0.0000 

0.0000 

-0.1367 

-0.1996 

-0.1273 

9 (degrees) 

0.0000 

0.0000 

2.6343 

3.8330 

3.8330 

density p 

1.6374 

0.3626 

2.1101 

1.8823 

0.4173 

pressure p 

0.3327 

0.3327 

0.4754 

0.4051 

0.4051 

sound speed c 

0.5333 

1.1333 

0.5616 

0.5489 

1.1658 

Mach number \M\ 

5.6250 

1.7647 

5.2956 

5.4396 

1.6335 


Table 4.3.1. Flow properties for the shock-wave/shear-layer test case in various regions of the 
flow: (1) upper stream inflow, (2) lower stream inflow, (3) upper stream after oblique 
shock, (4) upper stream after expansion fan, (5) lower stream after shock wave. 
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0.01, At = 0.01 K = 0.01, At = 0.02 



Fig. 4.1.4. Converting Vortex: Comparison of ACM66-ENT with the exact solution (solid line), 
illustrated by density profiles at the centerline y = 0, at I = 100,200,300 (curves 2 
- 4) for A t = 0.01,0.02,0.04 and k = 0.01, and for A t = 0.04 and k = 0.04 on a 
80 x 79 grid. 
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k — 0.01 on a 80 x 79 grid. The left most solid curve is the exact solution and the rest 
of the curves shifted to the right are the corresponding time sequences in an increasing 
order. 
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Density Contours (0.47 - 1 .04, 0.03 inc) 

CEN66, dt=0.02 
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Density Contours (0.47 -1 .04, 0.03 inc) 
CEN66, dt=0.01 
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Density Contours f0.47 -1 .04. 0.03 Inc) 
CEN66-ENT, dt=0.02 
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Fig. 4.1.10, Convecting Vortex: Comparison of CEN66-ENT with the exact solution (I.C.), 
illustrated by density contours at t = 50,100,200,300,370 for A t = 0.02 on a 
30 x 79 grid. 









Density Contours (0.47 - 1 .04, 0.03 inc) 

CEN66-ENT, dt=0.01 
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Fig. 4.1.1 1. Converting Vortex: Comparison of CEN66-ENT with the exact solution (I.C.), 
illustrated by density contours at t - 100,200,300,400,500 for A* = 0.01 on a 
80 x 79 grid. 








Density Contours (0.47 - 1 .04, 0.03 inc) 

ACM66, k=0.06, dt=0.04 



Fig. 4.1.12. Convccting Vortex: Comparison of ACM66 with the exact solution (I.C.), illustrated 
by density contours at t — 100,200,300,400,600, 600,700,800,000,1000,1100 for 
At = 0.04 and * = 0.06 on a 80 x 79 grid. 










Density Contours (0.47 - 1 .04, 0.03 inc) 
ACM66, k=0.06, dt=0.04 
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Density Contours (0.47 - 1 .04, 0.03 inc) 

ACM66-ENT, k=0.04, dt=0.04 











Density Contours (0.47 - 1 .04, 0.03 inc) 

ACM66-ENT, k=0.04, dt=0.04 
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Density Contours (0.47 - 1 .04, 0.03 inc) 
ACM66-ENT, k=0.01, dt=0.02 
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Density Contours (0.47 - 1 .04, 0.03 inc) 

ACM66-ENT, k=0.01 , dt=0.02 
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Fig. 4.1.14. (Cont.) 








Density Contours (0.47 - 1.04, 0.03 inc) 
ACM66-ENT, k-0.01, dt=0.01 
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Density Contours (0.47 - 1 .04, 0.03 inc) 

ACM66-ENT, k=0.01 , dt=Q.01 
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Density Contours (0.47 - 1 .04, 0.03 inc) 

ACM66S-ENT, k=0.005, dt=0.01 
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Fig. 4.1.16. Convecting Vortex: Comparison of ACM66S-ENT with the exact solu- 

tion (I.C.), illustrated by density contours at 1 = 100,200,300,400,500. 
600, 700, 800, 000, 1000, 1100 for A t = 0.01 and k = 0.005 on a 80 x 70 grid. 









Density Contours (0.47 - 1 .04, 0.03 inc) 

ACM66S-ENT, k=0.005, dt=0.01 
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Fig. 4.1.16. (Cont.) 
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Normalized Temperature Contours 

ACM66 & ACM66-ENT, k=0.7, k=0.35 for lin. field 

Uh-split Alpha=-5 Alpha=- 
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Fig. 4.2.3. Vortex pairing: Comparison of normalized temperature contours for a = -5 and -3 
with the un-split approach at time t = 160 on a 101 x 101 grid with k = 0.7 for the 
nonlinear fields and * = 0.35 for the linear fields using ACM66. 



Normalized Temperature Contours 

ACM66 & ACM66-ENT, k=0.7, none lin. field 



Fig 4 2.4. Vortex pairing: Comparison of normalized temperature contours for a = -3 with the 
un-split approach at time t = 180 on a 101 x 101 grid with i* = 0.7 for the nonlinear 
fields and * = 0 for the linear fields using ACM66. 
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Fig. 4.3.1. Schematic of the shock Impingement on a spatially-developing mixing layer. 


at t=120 (min= 0.314432 max 



Fig. 4.3.2 Shock-shear-layer interaction: The reference solution at t = 120 using ACM44. 

Contours are shown of (a) density, (b) pressure and (c) normalized temperature for a 
641 x 161 grid with * = 0.7 for the nonlinear fields and * = 0.35 for the linear fields. 






Shock Impingement on Mixing Layer 

Presure Contours (0.2 - 0.72, 0.02 inc) 



Fig. 4,3 3. Shock-shear-layer interaction: Comparison of the pressure contours for a = +10 with 
the un-split approach using ACM66 at t = 120 on a 321 x 81 grid with * = 0.7 for 
the nonlinear fields and k = 0.35 for the linear fields. 
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